Journal homepage    http://revistas.unitru.edu.pe/index.php/SSMM

Método de Diferencias Finitas Explícito para la solución de la ecuación de la onda acústica en un dominio heterogéneo.


Explicit Finite Differences Method for the solution of the acoustic wave equation in a heterogeneous domain.

Frank Henry Acasiete Quispe

Received, Feb. 04, 2019 , Accepted, Oct. 14, 2019 -

How to cite this article:
Acasiete, F. Método de Diferencias Finitas Explícito para la solución de la ecuación de la onda acústica en un dominio heterogéneo. Selecciones Matemáticas. 2019; 6(2):289-296. http://dx.doi.org/10.17268/sel.mat.2019.02.15

Resumen

Los Métodos de Diferencias Finitas (MDF) explícitos son usados para simulaciones numéricas en problemas de sísmica. En estos casos el medio es heterogéneo y para la garantía de una solución precisa y estable hay la necesidad de utilizar discretizaciones refinadas para las variables espacial y temporal, dando como resultado problemas de grandes dimensiones. En este trabajo presentamos el MDF explícito para la ecuación de onda en un dominio heterogéneo y mostramos los resultados a través de simulaciones computacionales.

Palabras clave. Ecuación diferencial parcial, onda, diferencias finitas, dominio heterogéneo.

Abstract

Explicit Finite Differences Methods (FDM) are used for numerical simulations in seismic problems. In these cases the medium is heterogeneous and for the guarantee of a precise and stable solution there is the need to use refined discretizations for spatial and temporal variables, resulting in large-scale problems. In this paper we present the explicit FDM for the wave equation in a heterogeneous domain and show the results through computer simulations.

Keywords. Partial differential equations, wave, finite differences, heterogeneous domain.

Introducción

La simulación de la propagación de ondas sísmicas, tanto en el Modelamiento Sísmico como en la Migración Reversa en el Tiempo, permite obtener información importante respecto a las estructuras geológicas que se encuentran en la sub-superficie. La solución numérica de los modelos matemáticos utilizados para la realización de estas simulaciones envuelve millones de grados de libertad, estos son determinados a cada instante de tiempo y para realizar estas simulaciones de forma eficiente se busca modelos numéricos sufucientemente robustos que sean estables y precisos, con un procesamiento computacional eficaz, que sea capaz de hacer viable los cálculos que son realizados

Los Métodos de Diferencias Finitas (MDF) explícitos son frecuentemente utilizados en la solución numérica de estos problemas y en este caso, el refinamiento de las discretizaciones espacial y temporal están relacionados a la menor longitud de onda que debe ser capturada y limitada por la condición de estabilidad del método utilizado, respectivamente [1,3,2]. La precisión de estos métodos depende de la orden de aproximación de los operadores de diferencias finitas y del refinamiento de la malla. Mallas más refinadas resultan en aumento de precisión, esto puede acarrear el crecimiento sustancial en el uso de memoria y tiempo de procesamiento, esto puede hacer que el costo computacional sea muy alto para el modelamiento de problemas realistas. En el caso de aplicaciones de MDF en dominios heteogéneos, como ocurre en problemas de sísmica, los algoritmos convencionales usan, para garantizar la estabilidad numérica de la solución, la discretización temporal dada por el menor paso de tiempo, definido por la sub-región del dominio con la mayor velocidad de propagación de la onda. De esta forma, existen diversos estudios en busca de una solución precisa y viable computacionalmente para este problema, entre ellas, esquemas de diferencias finitas de malla variable [9], utilización de operadores de alto orden para aproximar las derivadas espaciales y temporales [7,4] y por lo tanto la mayoria de los algoritmos utilizados para la resolución de este problema utilizan métodos explícitos, también existen trabajos que analizan la utilización de métodos implícitos [10,5].

La Ecuación de la Onda Acústica y el Método de Diferencias finitas

La ecuación de la onda acústica para problemas bidimensionales (2D), puede ser escrita como:

$$\frac{1}{c^{2}}\frac{\partial ^{2}u}{\partial t^{2}}-\left(\frac{... ...^{2}u}{\partial x^{2}}+\frac{\partial ^{2}u}{\partial y^{2}} \right) = f(x,y,t)$$ (1)

donde $u=u(x,y,t)$ es la amplitud de la onda en función de las variables espaciales en las direcciones $x$ y $y$, $t$ es la variable temporal, $c$ es la velocidad de propagación de la onda en el medio y $f(x,y,t)$ es el término fuente. Además de eso se dan las codiciones iniciales $u(x,y,0)=g(x,y)$ y $\frac{\partial u}{\partial t}(x,y,0)=h(x,y)$, $\forall x,y \in \Omega$ y condiciones de contorno adecuadas.

En este trabajo mostraremos la discretizacón del Método de Diferencias Finitas (MDF) expícito, con aproximaciones de segundo orden en el tiempo, de segundo orden y de cuarto orden para aproximar las derivadas espaciales. La aproximación explícita de diferencias finitas para la ecuación de la onda acústica utilizando aproximaciones de segundo orden en el tiempo y en el espacio, MDF(2-2), es dada por:

$$\frac{1}{c^{2}} \left( \frac{u_{i,j}^{t+1}-2_{u,j}^{t}+u_{i,j}^{t... ...+\frac{u_{i,j+1}^{t}-2_{u,j}^{t}+u_{i,j-1}^{t}}{\Delta y^{2}} \right)= f(x,y,t)$$ (2)

Desarrollando (2.2) obtenemos:

$$u^{t+1}_{i,j}=2u^{t}_{i,j}+C^{2}(u^{t}_{i+1,j}-2u^{t}_{i,j}+u^{t}_{i-1,j}-2u^{t}_{i,j}+u^{t}_{i,j-1}+u^{t}_{i,j+1})-$$

$$-u^{t-1}_{i,j}+ (c^{2} \Delta t^{2}f(x,y,t))+ O(\Delta t^{2}, h^{2}).$$ (3)

siendo $C$ el número de Courant $\left(C=\frac{c \Delta t}{h} \right)$, donde $\Delta t$ y $h$ definen las discretizaciones temporales y espaciales, respectivamente, con $\Delta x=\Delta y = h$. Los índices $i$ y $j$ designan los puntos de la malla espacial en las direcciones $x$ y $y$ respectivamente, en cuanto $t$ designa el instante del tiempo actual. Análogamente, se presenta en (2.5) la expresión del MDF explícito cuando se utilizan aproximaciones de diferencias finitas de segundo orden para la variable temporal y de cuarto orden para las variables espaciales, MDF(2-4).

$$\frac{1}{c^{2}}( \frac{u_{i,j}^{t+1}-2_{u,j}^{t}+u_{i,j}^{t-1}}{\... ...}+16_{i+1,j}^{t}-30u_{i,j}^{t}+16u_{i-1,j}^{t}-u_{i-2,j}^{t}}{12(\Delta y)^{2}}$$

$$- \frac{-u_{i,j+2}^{t}+16_{i,j+1}^{t}-30u_{i,j}^{t}+16u_{i,j-1}^{t}-u_{i,j-2}^{t}}{12(\Delta y)^{2}} ) = f(x,y,t)$$ (4)

Desarrollando (2.4) obtenemos:

$$u^{t+1}_{i,j}=2u^{t}_{i,j}+\frac{C^{2}}{12}[-u^{t}_{i+2,j}-u^{t}_{i-2,j}-u^{t}_{i,j+2}-u^{t}_{i,j-2} +$$

$$+16(u^{t}_{i+1,j}+u^{t}_{i-1,j}+u^{t}_{i,j-1}-60u^{t}_{i,j})]-$$

$$-u^{t-1}_{i,j} +(c^{2} \Delta t^{2}f(x,y,t)) + O(\Delta t^{2},h^{4}).$$ (5)

Las aproximaciones de Diferencias Finitas descritas anteriormente son condicionalmente estables y sus límites de estabilidad para dominios uni y bidimensionales, dados en términos del número de Courant, son presentados en la Tabla 2.1 que se encuentra en [8], para aproximaciones de segundo orden en el tiempo y de segundo y cuarto orden en el espacio.

Tabela 2.1: Límites de estabilidad

Dimensión	2do orden	4to orden
1D	$1$	$\frac{\sqrt{3}}{2}$
2D	$\frac{1}{\sqrt{2}}$	$\sqrt{\frac{3}{8}}$

Aplicaciones de estos algoritmos convencionales en medios heterogéneos, tiene la discretización espacial definida por la menor longitud de onda, que se desea representar y la discretización temporal constante, limitada por la condición de estabilidad de la sub-región que lleva el menor paso del tiempo, implicando en un gran costo computacional.

El error de dispersión

En el MDF los análisis de error de dispersión son realizados considerando un medio homogéneo y discretizaciones temporales y espaciales constantes. Explicaremos la influencia que esta discretización puede tener en dominios heterogéneos. Para esto, utilizamos la ecuación de la onda acústica unidimensional (1D).

$$\frac{1}{c^{2}} \frac{\partial ^{2} \varphi }{\partial t^{2}}= \frac{\partial ^{2} \varphi }{\partial x^{2}}$$ (6)

con condiciones iniciales y de contorno correctamente definidas.

Para el MDF explícito con aproximaciones temporal y espacial de segundo orden (2-2), se obtiene la siguiente expresión:

$$u^{t+1}_{i}=2u^{t}_{i}+C^{2}(u^{t}_{i+1}-2u^{t}_{i}+u^{t}_{i-1})-u^{t-1}_{i} + O(\Delta t^{2},h^{2})$$ (7)

y el error de dispersión, conforme presentado en [6], puede ser dado por:

$$\omega = \frac{2}{\Delta t}sen\left[ c \frac{\Delta t}{\Delta x}sen \left(\frac{\Delta x}{2} \right) \right] .$$ (8)

Colocando este resultado en términos de la velocidad de grupo normalizada, definida como $\frac{v_{g}}{c}=\frac{1}{c} \frac{d \omega}{dk}$, (2.8) puede ser reescrito como:

$$\frac{v_{g}}{c}= cos\left( \frac{k \Delta x}{2} \right) \left[ 1-... ...lta x^{2}} sen^{2} \left( \frac{k \Delta x}{2} \right) \right] ^{- \frac{1}{2}}$$ (9)

Para ejemplificar el comportamiento del error de dispersión cuando se utilizan las mismas discretizaciones, temporal y espacial, el número de Courant será alterado variándose la velocidad de propagación de onda.

Se puede observar que, como es presentado en [1], el error de dispersión en el MDF es menor cuando el número de Courant es mayor, en este caso próximo al límite de estabilidad del método.

Al haber una discontinuidad, como es el caso de nuestro problema, procederemos con nuestro método de diferencias finitas manteniendo nuestro número de Courant abajo del límite de estabilidad como está en la Tabla 2.1 para obtener los resultados que deseamos. A diferencia del método usado en [2] usaremos el MDF para nuestro problema y procuraremos obtener mejores resultados refinando más la malla del dominio usado.

Resultados Numéricos

Para evaluar los resultados obtenidos usando el método explícito, (2.1) es resuelto en un dominio heterogéneo 2D de $6000m$X$6000m$, subdividido en 2 sub-regiones, la división se encuentra en $x=3250m$.

Para nuestros resultados usaremos en los casos $c_{1}=V_{1},c_{2}=V_{2}=constantes$; la Sub-región 1 tiene la menor velocidad de propagación $V_{1}=1500m/s$, la Sub-región 2 tiene la mayor velocidad de propagación $V_{2}=2V_{1}=3000m/s$, con la siguiente condición inicial:

$$u(x,y,0)=g(x,y)=\frac{1}{12800 \pi} e^{ \left( - \frac{ \left( x^{2}+y^{2} \right)}{12800} \right)}$$ (10)

La Fig. 3.1 representa $g(x,y)$, cuya solución en $y=3000m$ es la siguiente:

Figura 3.1: Condición Inicial (Eq 3.1)

y $\frac{\partial u}{\partial t}=h(x,y)=0, \forall x,y$. La condición de contorno es la de Dirichlet, $u(\Gamma,t)=0,t>0$, siendo $\Gamma$ el contorno del dominio considerado, definido por $x \in [0;6000]$ y $y \in [0;6000]$.

Para las pruebas usamos fuente nula esto es $f(x,y,t)=0$. Usamos $\Delta t=0.008,0.004,0.002s$; $\Delta x =\Delta y = h=40,20,10m$, para mantener el número de Courant con el mismo valor y no tener problemas de dispersión en las simulaciones; usamos el método explícito de cuarto orden en el espacio y segundo orden en el tiempo. Los resultados son los siguientes; usamos las condiciones de frontera de Dirichlet, usamos estas condiciones para poder simular el caso en el cual el dominio no tuviera condiciones de frontera, ya que estamos en un dominio cuadrangular; para que la onda se propague uniformemente y el dominio tendrá reflexión, como se podrá ver en la Fig. 3.7:

Figura 3.2: Corte en $y=3000m$, para $t=0.6s$

Figura 3.3: Corte en $y=3000m$, para $t=0.9s$

Figura 3.4: Corte en $y=3000m$, para $t=1.2s$

En las Fig 3.2, 3.3, 3.4, tenemos el corte en la mitad del dominio en el eje $y=3000m$, con las diferentes discretizaciones espaciales y los tamaños de paso $\Delta t$.

Figura 3.5: Propagación de onda para $t=0.6s$

Figura 3.6: Propagación de onda para $t=0.9s$

Figura 3.7: Propagación de onda para $t=1.2s$

En las Fig. 3.5, 3.6, 3.7; tenemos la propagación de onda a través del dominio establecido y podemos observar la mayor propagación en donde la velocidad es mayor.

Conclusión

• La ecuación de onda es frecuentemente utilizada para modelar problemas en sísmica, que envuelven medios heterogéneos y el método de diferencias finitas explícito es ampliamente utilizado para encontrar su solución. Por ser un método condicionalemente estable se utiliza un incremento de tiempo constante en todo el dominio, dado por la situación más desfavorable, lo que implica en un gran esfuerzo computacional.

• En las Figuras 3.2, 3.3, 3.4; se muestra los cortes de la onda a través de instantes de tiempo en donde las curvas son más suaves a medida que se refina más la malla. En las Figuras 3.5, 3.6, 3.7; se muestra la propagación de onda, en donde se puede apreciar que ésta es más notoria en la zona con mayor velocidad de propagación.

• Otra variante para la resolución de estos métodos que permitan tener un ahorro considerable en el costo computacional y esto beneficie en la obtención rápida de resultados, a pesar de que su implementación requiera más esfuerzo y tiempo, serían los métodos adaptativos. Otra aplicación de los métodos Explícitos de Diferencias Finitas es para el caso en donde la ecuación (2.1) tenga término fuente esto es $f(x,y,t) \neq 0$; también se puede realizar simulaciones para dominios más complejos.

ORCID and License

Frank Acasiete Quispe https://orcid.org/ 0000-0002-3380-3606.

This work is licensed under the Creative Commons Attribution-NoComercial-ShareAlike 4.0.

Bibliografia

1

Alford, R.M., Kelly, K.R. and Boore, D.M. Accuracy of Finite-Difference Modeling of the Acoustic Wave Equation. Geophysics. 1974; 39:834-842.

2

Bartolo, L., Dors, C., and Mansur, W.J. A New Family of Finite Difference Schemes to Solve the Heterogeneous Acoustic Wave Equation. Geophysics. 2012; 77:187-199.

3

Bulcão, A. Modelagem e Migração Reversa no Tempo empregando operadores elásticos e acústicos. PhD thesis. 2004; COPPE-UFRJ, Brazil.

4

Chen, J. High-order time discretization in seismic modeling. Geophysics. 2007; 72:115-122.

5

Chua, C.,and Stoffab, P.L. Nonuniform grid implicit spatial finite difference method for acoustic wave modeling in tilted transversely isotropic media. Journal of Applied Geophysics. 2012; 76:44-49.

6

Cohen, G.C. Higher Order Numerical Methods for Transient Wave Equations. Springer. 2002.

7

Fornberg, B. The pseudospectral metrod-comparisons with finite differences for the elastic wave equation. Geophysics. 1987; 52:483-501.

8

Lines, L.R., Slawinski R., and Bording, R.P. A recipe for stability analysis of finite-difference wave equations computations. CREWES Research Report. 1998; 10:1-7.

9

Liu, Y., and Sen, M. Finite-difference modelling with adaptive variable-length spatial operator. Geophysics. 2011; 76:79-89.

10

Liu, Y., and Sen, M. A practical implicit finite-difference method: Examples from seismic modeling. Journal of Geophysics and Engineering. 2009; 6:231-249.

Notas al pie

... Quispe

Laboratorio Nacional de Computación Científica, MCTIC, Petrópolis- Brasil. (frankhaq@lncc.br).