Wave propagation: a finite difference modelling in a 3D fluid-solid configuration

Xiomara Contreras¹, Milagrosa Aldana²

¹ Departamento de Computación y TI. ² Departamento de Ciencias de la Tierra. Universidad Simón Bolívar. Valle de Sartenejas, Baruta, estado Miranda. xiomara@ldc.usb.ve, maldana@usb.ve

Abstract

In this work, the behavior of a wave field in a 3D solid-fluid configuration has been modeled. The elastodynamic equations that describe the problem were expressed, independently, in terms of both displacement- stress (EDE) and velocity-stress (EVE). Both systems were solved using a staggered-grid Finite Difference approach (SFD). A horizontal and an inclined surface were considered. In the first case, a transitional zone was used to solve the discontinuity at the interface. The other case introduces the medium heterogeneity into the grid, point by point, by means of the Lamé parameters and the density values. Our results indicate, for the first strategy, an increase in the amplitudes of the refractions and reflections; however, it does not modify the correct understanding of the geological model. The second strategy does not aggregate numerical effects to the results, but it could be unsuccessful in the presence of an irregular interface. Nevertheless, both strategies can be appropriate to solve basic 3D fluid-solid interface problems.

Keywords: propagation, waves, fluid-solid, finite-differences, transition-zone.

Propagación de ondas: modelado por diferencias finitas en una configuración 3D fluido-sólido

Resumen

En el presente trabajo se modela el comportamiento de un campo de ondas en una configuración tridimensional fluido-sólido mediante las ecuaciones elastodinámicas expresadas en desplazamiento-esfuerzo (EDE) y en velocidad-esfuerzo (EVE), respectivamente. Ambos, EVE y EDE, se resuelven en forma aproximada por diferencias finitas sobre una staggered-grid (SDF). La simulación numérica se aborda con ambos métodos SDF EVE y EDE en dos medios con contacto fluido-sólido, donde el primero presenta una interfase horizontal y el otro una interfase inclinada. En el primer caso la interacción entre el fluido y el sólido se establece mediante una zona de transición posicionada entre ambos. En el otro caso se introduce la heterogeneidad en el medio, punto a punto en la malla numérica, mediante los parámetros de Lamé y los valores de la densidad. Los resultados obtenidos muestran que la refracción y reflexión del campo de ondas parece incrementada al considerar la zona de transición pero sin modificar la correcta interpretación del modelo geológico. Con el otro tratamiento no se observan efectos atribuibles al modelado numérico de la interfase, pero la estrategia puede no ser efectiva en presencia de una interfase irregular. Ambas estrategias pueden ser apropiadas al resolver el problema indicado en interfases fluido-sólido, básicas, en 3D.

Palabras clave: propagación, ondas, fluido-sólido, diferencias-finitas, zona-transición.

Recibido el 29 de Enero de 2011

En forma revisada el 5 de Marzo de 2012

Introducción

Los métodos de exploración sísmica, ampliamente utilizados en la caracterización de yacimientos, tienen como base el estudio de la propagación de ondas acústicas generadas artificialmente [1]. En su desplazamiento en el subsuelo, las ondas atraviesan capas con propiedades físicas distintas, pudiendo interactuar con interfases sólido-sólido o fluido-sólido. Si bien las primeras han sido tratadas ampliamente, la complejidad del segundo caso ha limitado su estudio, principalmente en 3D. Sin embargo, el interés de la exploración sísmica en interfases fluido-sólido ha aumentado con la exploración sísmica marina [2]. De ahí la importancia de estudiar el tratamiento de este tipo de interfases en medios 3D. En general, una interfase se modela matemáticamente mediante las condiciones explícitas de borde en ésta, pero tales condiciones pueden ser numéricamente complejas de establecer para interfases con una geometría irregular. Sin embargo, existe la alternativa de modelar implícitamente los cambios presentes en la proximidad de la interfase [2]. En el presente trabajo se aplican, en 3D, dos formulaciones de las ecuaciones elastodinámicas para describir matemáticamente el comportamiento de un campo de ondas al desplazarse entre dos regiones con contacto plano, una acústica (fluido) y una elástica (sólido). Las ecuaciones se expresan tanto en términos de desplazamientos-esfuerzos (EDE) como de velocidades- esfuerzos (EVE). Las derivadas parciales, en los sistemas EDE y EVE, se determinan aproximadamente por diferencias finitas sobre una staggered-grid.

La interacción entre una interfase fluido- sólido y un campo de ondas ha sido investigada en 2D, entre otros, por van Voosen et al. [3]. Estos autores emplean un método staggered-grid en diferencias finitas EVE, y una estrategia explícita denominada “Images method” para identificar la respuesta de una interfase horizontal. Komatitsch et al. [4] investigan una interfase, horizontal y sinusoidal, respectivamente, mediante la conjunción de las condiciones de borde explícitas en la interfase con un método de elemento espectral SEM (spectral element method). En medios 3D, Käser y Dumbser [5] investigan la propagación de ondas en un medio fluido-sólido mediante un método de elementos finitos, tratando explícitamente una interfase sinusoidal discontinua en un medio marino. Contreras y Aldana [6], por su parte, aplican un método EVE staggered- grid en diferencias finitas y tres estrategias, dos implícitas y una explícita, para resolver la interacción ondas-interfase en un medio fluido-sólido con contacto plano y horizontal. El presente trabajo complementa y extiende el enfoque aplicado en el de Contreras y Aldana [6] utilizando, dos métodos y dos estrategias implícitas, en diferencias finitas, y dos modelos geológicos fluido- sólido con contacto plano. En uno de los modelos la interfase es horizontal y en el otro está inclinada. La interfase horizontal se modela aquí mediante una zona de transición constituida por múltiples capas sólidas para suavizar los cambios entre las dos regiones. En el caso de la interfase inclinada, se introducen los cambios estructurales punto a punto en la malla numérica utilizando los parámetros de Lamé del medio, por lo que no se realiza tratamiento alguno para suavizar el contraste existente.

Planteamiento y solución del problema

Los sistemas EDE y EVE, en 3D, se plantean mediante un sistema lineal de nueve ecuaciones, cada uno. El planteamiento para EVE aparece en Contreras y Aldana [6]. Se muestra aquí el correspondiente a EDE mediante tres ecuaciones representativas (las restantes son análogas) y las hipótesis siguientes: un volumen rectangular V, acotado, con frontera externa G que representa el medio de propagación; un subconjunto W de V, abierto y acotado, W = W_F È W_S donde, W_F y W_S son dominios abiertos en ambos fluido y sólido, respectivamente;

son espacios de funciones cuadrado integrables [7].

Para cada t, (0 < t < T) se determinan en H₀¹ (W) las funciones incógnitas u_z, u_x , u_y, s_zz, s_xx, s_yy, s_zx, s_zy, s_xy que satisfacen en  las ecuaciones (1) a la (3), además de las seis ecuaciones restantes que suponemos implícitas:

En estas ecuaciones, [0, T] es un intervalo real; t es la variable temporal (t Î [0, T] ); x es un vector en R³; 1/r es el grado de iluminación dado por el inverso de la densidad r; l y μ son los parámetros de Lamé; f_z, f_x, f_y: son las componentes de la fuente; u_x , u_y, u_z son las componentes del desplazamiento; s_zz, s_xx, s_yy, s_zx, s_zy, s_xy son los esfuerzos normales y deformantes; son las respectivas velocidades de ondas P y S.

Las condiciones iniciales y de borde en la interfase fluido-sólido son las mismas para ambos métodos EVE y EDE: el medio está en equilibrio en t=0, es decir desplazamientos y esfuerzos son cero en t=0; en la frontera , para cada t, las funciones solución se anulan; en la interfase fluido- sólido, los desplazamientos y los esfuerzos normales son continuos y los esfuerzos deformantes se anulan [3].

Los sistemas en EVE y EDE son resueltos independientemente mediante métodos, a cuarto orden de aproximación espacial y a segundo orden de aproximación temporal, en diferencias finitas sobre una staggered-grid. En ambos casos se emplea la misma definición para la staggered-grid, reemplazando apropiadamente los desplazamientos por las velocidades [6]. La aproximación por diferencias finitas para las ecuaciones (1) a la (3) es la siguiente:

Larsen [8] aplica, en simulaciones con modelos terrestres, la fórmula de estabilidad D_t < (0,494* D_z)/V, donde V es la velocidad máxima de la onda P en el medio. En Contreras y Aldana [9] fue analizada la condición de estabilidad para el método EVE en una configuración fluido-sólido. En ese trabajo, se varían el time-step y la escala en la malla (malla fina y malla gruesa). La variación del time-step confirmó la condición de estabilidad heterogénea con resultados estables, mientras que el cambio de escala espacial permitió identificar la presencia de artefactos numéricos originados por la representación de la fuente en una malla gruesa. Sin embargo, en una malla fina se identificó, en imágenes de alta calidad, la interacción del campo de ondas con una interfase fluido-sólido con escasos artefactos numéricos y observándose la reflexión-transmisión de la onda P y las ondas S convertidas en la zona de transición. La respuesta de la interfase fue acentuada de forma positiva por una zona de transición de poco espesor.

Simulación

Para la simulación se consideró un volumen rectangular, acotado y dividido por una interfase plana, horizontal o inclinada. Las condiciones de borde en la interfase se modelan implícitamente en cada caso. En el presente trabajo no se aplican las condiciones explícitas de borde en la interfase. En Contreras y Aldana [6] se describen esas condiciones en diferencias finitas para la interfase horizontal. La zona de transición multicapa se define de acuerdo a lo siguiente: en cada capa la densidad y las velocidades sísmicas se calculan mediante la fórmula lineal p_f+[i*(p_s-p_f)/N], 0 ≤ i ≤ N, donde N es el número de capas y p_f, p_s representan el valor del parámetro considerado (densidad, velocidades de onda P o S), en el fluido y en el sólido, respectivamente. Es importante destacar que en Contreras y Aldana [6], una de las estrategias implícitas que describe el contacto fluido-sólido, inserta una zona de transición en el medio de propagación, soldada y paralela a la interfase, y constituida por una capa sólida homogénea.

Como se indicó anteriormente, las simulaciones se realizan considerando las dos formulaciones, EVE y EDE, una interfase horizontal y otra inclinada y dos tipos de sólidos. El fluido se coloca en la parte superior del volumen; la malla numérica es uniforme y su dimensión corresponde a 200 x 250 x 200, para la interfase horizontal, y 200 x 200 x 200 para la interfase inclinada; la fuente se ubica en el centro del fluido y se dispara en los instantes D_t + (D_t /2) para EVE y D_t para EDE (D_t identifica el time-step). El espesor de la zona de transición es N (N=20D_z). Se utilizan bordes absorbentes sobre todas las caras del volumen [10]. Estos bordes se intersecan entre ellos a lo largo de los ejes y esquinas del volumen, atenuando o anulando las reflexiones del campo de ondas en dichas intersecciones. El tiempo total de cada simulación es 700D_t (3,451mseg). Se monitorea el transcurrir de la propagación usando la solución para s_zz en instantes de tiempo indicados, en cada caso. A continuación se presentan dos problemas de estudio. En cada caso se considera un sólido diferente, mientras que el fluido es el mismo. Las imágenes iniciales de la propagación se omiten en todos los casos.

Problema 1

Los valores de las velocidades y densidades se seleccionaron de van Voosen et al. [3]. En el fluido (un líquido [11]) la velocidad de onda P es 1500 m/seg y la densidad 1000 Kg/m³. En el sólido elástico [3], las velocidades de ondas P y S son 3500 m/seg y 2000 m/seg, respectivamente, y la densidad es de 2500 Kg/m³. El espaciamiento en los ejes es de 0,035m y 0,00493mseg es el time-step.

En las Figuras 1 y 2 se muestra, para la interfase horizontal, el resultado de la propagación en los primeros pasos en tiempo. Se observa una onda P esférica que se expande en el fluido hasta alcanzar la interfase fluido-sólido, a partir de la cual se expande un frente de onda en el sólido. Se identifica la reflexión de la onda P, al inicio de la zona de transición, sobre la interfase fluido-sólido, y en el borde superior de la malla con atenuación. En la Figura 3, las ondas P reflejadas en la zona de transición y sobre la interfase fluido-sólido se fusionan visualmente en un campo de ondas de amplio espesor que asciende hacia el borde superior. Se observa también la intersección entre la onda P reflejada en el borde superior y el campo de ondas visualmente fusionado. Por otra parte, la onda P transmitida al sólido toca el borde inferior del volumen, mientras que la onda S, desarrollada en el sólido, está alcanzando ese borde. En el caso de EDE, se observa un arco sombreado que asciende desde el borde inferior, el cual podría estar asociado a la onda P reflejada y atenuada allí por la acción del borde absorbente. En el caso de EVE, se observa también un arco de color blanquecino, ascendiendo desde el borde inferior. En la Figura 4, el frente de ondas alcanza el borde superior en ambos EDE y EVE. Se observa en EVE la atenuación de las ondas provenientes del fluido que aún son transmitidas a la zona de transición; la onda P reflejada en el borde inferior de la malla ha tocado la interfase y regresa a la zona de transición. Para EDE se observa algo similar intensificado por la presencia de una acumulación de luz bajo la interfase fluido-sólido; el grado de iluminación puede estar asociado al inverso de la densidad. Las imágenes para EDE presentan bordes suaves mientras que, para el método en EVE, son delgados y bien definidos. Esta diferencia parece tener relación con el desfase de D_t/2 entre los tiempos de disparo de la fuente para EVE y EDE.

Respecto a la interfase inclinada, en la Figura 5, el campo de ondas está contenido en el fluido. En la Figura 6, la onda P ha interactuado con la interfase identificando la inclinación de ésta; se observa, además, la refracción y la reflexión del frente de ondas allí. En la Figura 7, observamos dos frentes de ondas avanzando en el sólido; el más adelantado corresponde a la onda P transmitida y el que le sigue corresponde a la onda S convertida. En el fluido, por otra parte, observamos la intersección entre la onda P reflejada en la interfase y la onda P reflejada en el borde superior de la malla. En la Figura 8, en el fluido, se observa respecto a la onda P que la directa ha tocado el borde derecho, las respectivas reflexiones en la interfase y en el borde superior han tocado respectivamente el borde superior y la interfase. En el sólido, se observa la reflexión de la onda P en el borde inferior y con apariencia atenuada en el borde izquierdo; la onda S está tocando el borde inferior.

La Figura 9 muestra los sismogramas para ambos métodos. En los dos casos se observan dos hipérbolas centradas, separadas en el tiempo, asociadas a la onda P directa y a su reflexión en el borde superior de la malla; una semi hipérbola desplazada hacia la izquierda está asociada a la reflexión de la onda P en la interfase. Más a la izquierda se identifican semi arcos asociados a las reflexiones de las ondas P y S sobre el borde izquierdo de la malla.

Problema 2

Los valores en el fluido son los mismos utilizados en el problema 1. El sólido corresponde a una roca sedimentaria detrítica (una lutita). Una trampa de hidrocarburos se caracteriza por la presencia de una capa porosa y permeable, como una arenisca, donde se acumula el fluido, y una capa impermeable que actúa como sello e impide la migración del hidrocarburo. Las lutitas generalmente actúan como roca sello [12]. Este tipo de litología se encuentra en regiones inundadas y en deltas de los ríos, entre otros ambientes [13]. Por lo tanto, este caso es de interés en prospección y, en este trabajo, lo definimos y analizamos hipotéticamente. Los valores para las velocidades P, S y la densidad son 2900 m/seg, 1330 m/seg, y 2290 Kg/m³, respectivamente [14]. El espaciado en la malla es 0,029 m y se utiliza el mismo valor del time-step que en el problema 1. Respecto a la interfase horizontal, en la Figura 10 se observan efectos atenuados de reflexión que ocurren en la zona de transición. En esta zona se observa un frente de ondas que asciende hacia el borde superior, el cual parece la unificación visual asociada, entre otros, a la onda P reflejada en las capas de la zona de transición. Al mismo tiempo puede observarse la reflexión de la onda P en el borde superior y un grueso frente de ondas transmitido a la región ocupada por la lutita. Hay efectos de iluminación sobre este último campo de ondas en ambos EDE y EVE, ligeramente acentuado en el caso de EDE. En la Figura 11, la onda S está llegando al borde inferior mientras que la reflexión de la onda P allí está atenuada por la acción de los bordes absorbentes. En la Figura 12 los efectos que se observan, en la zona de transición, pueden atribuirse a la refracción atenuada, de la onda P previamente reflejada en el borde superior de la malla, al hacer contacto con las capas sólidas de la zona.

Respecto a la interfase inclinada, la expansión del campo de ondas se comporta de forma similar a lo observado en el problema 1. En la Figura 13 observamos en el sólido el aspecto nítido de la onda P transmitida, de la onda S convertida allí y de la reflexión de la onda P en el borde izquierdo. Esto difiere de lo observado en el sólido del problema 1 y puede atribuirse a una mayor atenuación de la señal en ese caso. En las Figuras 14 y 15 observamos la intersección entre diferentes campos de ondas, más acentuado en la última.

Conclusiones

Los métodos en diferencias finitas, EDE y EVE, utilizados en el presente estudio reportan información mutuamente complementaria en relación al comportamiento del campo de ondas al viajar entre los dos medios de interés, la reflexión y la transmisión de las ondas directas entre fluido y sólido, las ondas convertidas en las regiones sólidas y la identificación, en cuanto a posición e inclinación, de la interfase. La zona de transición, tal como fue implementada, agrega una reflexión y refracción extra del campo de ondas en la vecindad de la interfase horizontal, como pudo identificarse en ambos métodos. Sin embargo, este efecto numérico no parece desviar la correcta interpretación del modelo geológico. La estrategia aplicada para la interfase inclinada no reporta efectos numéricos, pero no disminuye los efectos que la interfase fluido-sólido introduce, lo que pudiera ser desventajoso en su aplicación en configuraciones similares con contrastes severos. Ambas estrategias resuelven adecuadamente la interacción ondas-interfase fluido-sólido en los casos de prueba considerados.

Agradecimientos

A los operadores del laboratorio LDC, USB, por mantener operativo el ambiente de cómputo que requirió la fase experimental de este trabajo.

Referencias bibliográficas

1. De Hoop A. T., Van der Hijden J. H. T.: “Generation of acoustic waves by an impulsive line source in a fluid / solid configuration with plane boundary”. J. Acoust. Soc. Am., Vol. 74, Nº 1 (1983) 333-342.        [ Links ]

2. Vireux J.: “P-SV wave propagation in heterogeneous media: Velocity-stress finite difference method”. Geophysics, Vol. 51 (1986) 889-901.        [ Links ]

3. Van Voosen R., Robertsson J. O. A., Chapman Ch. H.: “Finite-Difference Modeling of wave propagation in a fluid-solid configuration”. Geophysics, Vol. 67, (2002) 618-624.        [ Links ]

4. Komatisch D., Barnes C., Tromp J.: “Wave propagation near a fluid-solid interface: A Spectral element approach”. Geophysics, Vol. 65, Nº 2, (2000) 623-631.        [ Links ]

5. Käser M., Dumbser M. “A highly accurate discontinuous Galerkin method for complex interfaces between solids and moving fluids”. Geophysics, Vol. 73, Nº 3, (2008) T23-T35.        [ Links ]

6. Contreras X., Aldana M.: “Three strategies to solve the wave propagation problem in a 3D fluid-solid interface”. Revista de Ingeniería, UCV, Vol. 4, (2009).        [ Links ]

7. Johnson C.: “Numerical solution of partial differential equations by the Finite Element method”. Cambridge University Press. Sweden. 1994.        [ Links ]

8. Larsen Sh.:“E3D: 2D/3D finite-difference wave propagation code”. Livermore National Laboratory. Livermore, CA 94551. http://crack.seismo.unr.edu/ftp/pub/louie/class/455/e3d/e3d.txt. 1992.        [ Links ]

9. Contreras X. and Aldana M.: “Wave propagation in a 3D fluid-solid configuration: Staggered- grid finite difference modeling and stability analysis” Vol. 30, (2011) SM P1.        [ Links ]

10. Cerjan Ch., Kosloff D., Kosloff R., Reshef M.: “A nonreflecting boundary condition for discrete acoustic and elastic wave equations”. Geophysics, Vol. 50, (1985) 705-708.        [ Links ]

11. Levander A.: “Fourth-order Finite Difference P-SV seismograms”. Geophysics, Vol. 53, Nº 1, (1988) 1425-1436.        [ Links ]

12. Sheriff R. and Geldart Ll.:“Exploration Seismology”. Cambridge University Press, Second edition. United States of America. 1995.        [ Links ]

13. Auque L., Gimeno M., Gómez J.: “Rocas Siliciclásticas II. Lutitas”. Grupo de modelización Universidad de Zaragoza. 2010. http://gmg.unizar.es/gmgweb/Asignaturas/ExogenaII/Teoria/Tema_4.        [ Links ]

14. Upadhyay S. K.: “Seismic Reflection Processing: With special reference to anisotropy”. Geophysics & Geodesy. Springer. Germany. 2004.        [ Links ]