Este artículo muestra la resolución del problema de propagación de ondas sísmicas en 2-D, mediante la utilización de esquemas explícitos en diferencias finitas generalizadas (GFD), lo que permite la utilización de mallas regulares e irregulares.
Puesto que se utiliza un método explícito, es necesario obtener la condicióń de estabilidad, lo que se ha realizado mediante un análisis de von Neumann. También se obtienen las relaciones de dispersión en la estrella de las velocidades de fase para las ondas P y S, así como las de las velocidades de grupo.
Dada la importancia que en la aplicación del método tiene el control sobre la irregularidad de la malla, se han definido unos índices de irregularidad para la estrella (IIS) y la malla (IIC), analizándose su relación con la dispersión y el paso de tiempo utilizado en los cálculos.
This paper shows the solution to the problem of seismic wave propagation in 2-D using generalized finite difference (GFD) explicit schemes. Regular and irregular meshes can be used with this method.
As we are using an explicit method, it is necessary to obtain the stability condition by using the von Neumann analysis. We also obtained the star dispersion formulas for the phase velocities for the P and S waves, as well as the ones for the group velocities.
As the control over the irregularity in the mesh is very important in the application of this method, we have defined an index of irregularity for the star (IIS) and another for the cloud (IIC), analyzing its relationship with the dispersion and time step used in the calculations.
Son muchos los trabajos publicados sobre propagación de ondas sísmicas, incluso centrándonos en los relativos a la utilización del método de Diferencias Finitas (FD) y dentro de ellos los que tratan aspectos relativos al tipo de malla, la estabilidad y dispersión, que son los aspectos en que fundamentalmente se centra el presente trabajo. Así, durante más de diez an˜os se utilizarón mallas convencionales en FD, pudiéndose destacar los trabajos de Altermann y Karal [1] y Kelly et al. [2].
Los esquemas Staggered-grid (S-G) en FD fueron introducidos por Madariaga [3], popularizándose su uso en problemas de propagación de ondas con los trabajos de de Virieux [4], [5] en los que se usa un esquema de segundo orden para la formulación en velocidad-tensión con ondas SH y P-SV respectivamente.
Otros esquemas han sido presentados por Levander [6] para una formulación en cuarto orden para ondas P-SV, Graves [7], Graves et al. [8], Ohnimato and Chovet [9] que lo aplican en un esquema de segundo orden en la formulación desplazamiento-tensión. También pueden destacarse los trabajos de Kristek et al. [10], [11].
A pesar de los trabajos citados de Virieux y Levander sobre dispersión y Crase et al. [12] que obtienen condiciones de estabilidad para ondas P-SV para un orden arbitrario de aproximación, Moczo, Kristek y Halada [13] se dan cuenta de que la dispersión de los esquemas S-G en 3-D no están suficientemente investigados y lo estudian para esquemas de cuarto orden en el espacio y segundo orden en el tiempo con validez para formulaciones en desplazamiento-tensión, velocidad-tensión y desplazamiento-velocidad-tensión. También Graves et al. [8] han estudiado la estabilidad y Mozco et al. [14] presentan trabajos recopilatorios de gran interés.
Las mallas parcialmente Staggered(P-S-G) fueron utilizadas por primera vez por Andrews [15] para medios anisótropos y Zhang [16] las utiliza para la formulación velocidad-tensión en 2-D. Cruz-Atienza y Virieux [17] las utilizan para la simulación de una rotura de falla en 2-D.
El primero en utilizar una malla con espaciado variable fue Boore [18] para un caso veritcal 1-D, Mikuano y Miyatarke [19] usan un espaciado variable en un medio elástico en 2-D para analizar el proceso de rotura de falla. Pitarka [20] presenta un esquema de cuarto orden para la formulación en velocidad-tensión sobre una malla rectangular con espaciado variable y Wang et al. [21] combinan mallas con diferentes espaciados.
Los métodos que en general se pueden denominar métodos sin malla, surgieron hace aproximadamente treinta an˜os, aunque el interés y el esfuerzo investigador ha sido mínimo hasta hace poco tiempo. Así pues, la referencia más alejada en el tiempo y que puede condiderarse como punto de partida es el método de partículas (SPHM) que fue desarrollado por en 1977 por Lucy [22]. En 1982 y 1988 Monaghan ([23] y [24]) desarrolla una explicación más rigurosa del método.
Una camino paralelo al indicado en el desarrollo de aproximaciones sin malla se basa en la idea de realizar una aproximación local mediante el método de mínimos cuadrados móviles de Lancaster y Salkauskas (1981) ([25]). Nayroles, Touzot y Villon (1992) ([26]) fueron los primeros en utilizar una aproximación de este tipo en un método de Galerkin, denominado método de los elementos difusos (DEM).
Belytschko, Lu y Gu (1994) ([27]) refinaron y desarrollaron una implementación alternativa del método de los elementos difusos, que mejoró la precisión al utilizar multiplicadores de Lagrange para imponer las condiciones de contorno esenciales y un orden de cuadratura mayor en la integración. Ellos denominaron al método Galerkin libre de elementos (EFG).
Un salto importante en el desarrollo de estos métodos se ha producido a partir de los trabajos de Duarte y Oden (1995) ([28]) y de Babuska y Melenk (1995) ([29]) que realizan una generalización muy interesante de la aproximación por mínimos cuadrados móviles usando el concepto de partición de la unidad, que los primeros denominaron nubes hp (hp clouds) y de partición de la unidad (PUFEM) los segundos.
Otro camino en la evolución de los métodos sin malla se ha constituido mediante el desarrollo del método de diferencias finitas generalizadas (GFDM) que surge como consecuencia de la evolución del método de diferencias finitas clásicas. Los trabajos de Jensen (1972) ([30]) y Perrone y Kao (1975) ([31]) pusieron las bases del método. Liszka y Orkisz (1980) ([32]) y Orkisz (1998) ([33]) han realizado importantes contribuciones a la mejora y desarrollo del método.
El método de diferencias finitas generalizadas (GFDM) se basa en la utilización de mínimos cuadrados móviles, lo que permite aplicar esquemas en diferencias finitas en dominios irregulares. Al desarrollo de este método han contribuido los trabajos de los autores. Así, se han obtenido fórmulas explícitas en diferencias para el caso de estrellas irregulares, analizandose la influencia de los parámetros fundamentales que aperecen en la formulación [34], se ha obtenido un método adaptativo para la mejora de la solución aproximada utilizando GFDM [35], en [36] se muestra la aplicación del método adaptativo en GFDM a la resolución de ecuaciones diferenciales en derivadas parciales de segundo orden y [37] es un resumen global de GFDM que incluye las aportaciones realizadas por los autores y su aplicación a problemas elásticos.
Los artículos [38] y [39] muestran la aplicación del método de diferencias finitas generalizadas a la resolución de ecuaciones en derivadas parciales dependientes del tiempo.
Los métodos numéricos constituyen una herramienta irremplazable en la investigacíón de los movimientos sísmicos y entre ellos el actualmente dominante es el DF, es más Moczo et al. [40] opinan que está por llegar el mejor momento de la aplicación del método Df en sismología. Tras estas consideraciones parece oportuno presentar las aplicación del GFDM a este tipo de problemas.
Con el presente artículo se pretende establecer una base sólida para la aplicación de un método que ofrece grandes posibilidades entre las que inicialmente hay que destacar la capacidad de utilizar mallas irregulares directamente, evitando artificios en la realización de algunos modelos híbridos. Se desea hacer hincapié en la solidez metodológica puesto que no solo se presenta la formulación en GFD para la aplicación indicada, sino que se estudia de forma rigurosa su estabilidad y se analiza detalladamente la dispersión.
2Esquema en diferencias finitas generalizadas para el problema de propagación de ondas sísmicas en dominios homogéneos2.1Planteamiento del problemaEl problema de propagación de ondas sísmicas en 2-D viene caracterizado por el sistema de ecuaciones diferenciales en derivadas parciales
con las condiciones inicialesy la condición de contornosiendo f1(x, y), f2(x, y), f3(x, y), f4(x, y), g1(t) y g2(t) funciones conocidas,ρ es la densidad, λ cortante de Lamé y μ el módulo de rigidez. Las funciones U(x, y, t), V(x, y, t) son al menos dos veces diferenciables con continuidad en el dominio Ω con frontera Γ.2.2Esquema en diferencias finitas generalizadas
Para la obtención de las fórmulas explícitas en diferencias finitas de las derivadas espaciales, una vez discretizado el dominio Ω∪Γ, a cada nodo del mismo se le asigna un conjunto de nodos a su alrededor que se denomina estrella. En la figura 1 se muestra una estrella irregular con ocho nodos seleccionados con el criterio de los cuatro cuadrantes (dos nodos por cuadrante, que son los más próximos al nodo central, considerando este como origen).
La derivada segunda respecto del tiempo en el nodo central de la estrella se aproxima por:
donde u0n, v0n, u0n+1 y v0n+1 son los valores aproximados de las funciones U(x, y, t), V(x, y, t) en el nodo central de coordenadas espaciales (x0, y0) para los tiempos n▵t y (n+1)▵t respectivamente.
Las derivadas espaciales se aproximan por las expresiones explícitas en diferencias finitas generalizadas, obtenidas en [36] y [38]:
donde N es el número de nodos de la estrella cuyo nodo central tiene de coordenadas (x0, y0) (en este trabajo N=8 y son seleccionados mediante el criterio de los cuatro cuadrantes), m0, η0 y ζ0 son los coeficientes que multiplican a los valores aproximados de las funciones U y V en el nodo central para el tiempo n▵t (u0n y v0n respectivamente) de las expresiones explícitas en diferencias finitas generalizadas para las derivadas espaciales ∂2∂x2, ∂2∂y2 y ∂2∂x∂y respectivamente, mj, ηj y ζj son los coeficientes que multiplican a los valores aproximados de las funciones U y V en el resto de los nodos de la estrella para el tiempo n▵t (ujn y vjn respectivamente) de las expresiones explícitas en diferencias finitas generalizadas para las derivadas espaciales ∂2∂x2, ∂2∂y2 y ∂2∂x∂y respectivamente ([32], [33], [34], [37]).
Si en el sistema de ecuaciones 1 se sustituyen las derivadas parciales por las expresiones 4 y 5 se obtiene:
que constituyen el esquema en diferencias finitas generalizadas para el pro-blema planteado.3Criterio de estabilidad
Para estudiar la estabilidad se utiliza el análisis de von Neumann. Para lo cual se realiza una descomposición armónica de la solución aproximada de la forma
donde x0 es vector de posición del nodo central de la estrella, xj, j=1, …, N son los vectores de posición del resto de los nodos de la estrella y hj son los vectores de las posiciones relativas de los nodos de la estrella con respecto al nodo central cuyas coordenadas son hjx=xj−x0, hjy=yj−y0.
ξ=e−iw▵t es el llamado factor de amplificación, cuyo valor determinará la condición de estabilidad, w es la frecuencia angular y k es el vector del número de onda (ver figura 2)
Al sustituir las expresiones 7 en la ecuación 6, y simplificar, se obtiene.
Teniendo en cuenta que los coeficientes mj, ηj, ζj, son particiones de la unidad, se tieneal sustituir las expresiones 9 en la ecuación 8, se obtiene el sistema de dos ecuaciones con dos incógnitas, A y B, siguienteen el que al despejar B de la segunda ecuación del sistema 10, y sustituir en la primera ecuación, y tras operar, se obtiene
Tras operar, se obtienen las siguientes condiciones:
Parte real
Parte imaginaria
Al operar con estas ecuaciones, simplificando mediante el uso de criterios conservadores y eligiendo entre ambas la condición más restrictiva, se obtiene
que es la condición de estabilidad de la estrella. Por tanto, el paso de tiempo (▵t) deberá ser el mínimo de los obtenidos para cada una de las estrellas de la malla.4Dispersión en la estrella4.1Dispersión de las ondas P
De la ecuación 12 se tiene
dondedondeySe sabe quedonde cgrid y λgrid son la velocidad de fase (αgrid o βgrid) y la longitud de onda (λPgrid y λSgrid) en la estrella respectivamente ([4], [5], [6], [7], [12], [13], [14]).
Si se definen las relaciones:
La razón de velocidades esy por tantoSi se sutituyen las ecuaciones 15, 20, 21 y 23 en la ecuación 18 particularizada para las ondas P, tras operar, se obtiene la relación de dispersión para las ondas P4.2Dispersión de las ondas S
Se sabe que
Si se sustituyen las ecuaciones 15, 19 y 21 en la ecuación 25, tras operar, se obtiene4.3Dispersión de la velocidad de grupo ondas P
Por definición la velocidad de grupo para las ondas P, se obtiene al derivar w, dada por la expresión 15, respecto de k, por tanto
dondedesignando4.4Dispersión de la velocidad de grupo ondas S4.5Irregularidad de la estrella y dispersión
Los coeficientes, m0, η0, ζ0, que aparecen en el criterio de estabilidad dependen del número de nodos en la estrella, de las coordenadas de cada nodo respecto del nodo central de la estrella y de la función de ponderación elegida ([34]).
Fijados el número de nodos por estrella, en nuestro caso 9 (N=8), y la función de ponderación, Ω(hjx,hjy)=1(hjx2+hjy2)3, la expresión
depende de la razón de velocidades y de las coordenadas relativas de cada nodo de la estrella respecto del nodo central de la misma.
Al disminuir las coordenadas relativas de los nodos de la estrella, aumentan los coeficientes m0, η0, ζ0, lo cual ocurre al aumentar el número de nodos en la malla. Para un mismo número de nodos en la malla el mínimo valor de los coeficientes mencionados se obtiene para la malla regular ([36]).
Sea τl la media de las distancias de los nodos de la estrella l a su nodo central. Sea τ la media de todas las τl de las estrellas de la malla, una vez discretizado nuestro dominio. Se tiene que:
en las que la barra superior indica valores medios.
La condición de estabilidad puede reescribirse
En el caso de malla regular de ocho nodos la expresión 35 resulta serSi se multiplica el segundo miembro de la expresión 36 por el factorse obtiene la expresión 35.
El factor obtenido en 37, que para una estrella regular (esquema de ocho nodos más el nodo central) es la unidad, resulta ser un indicador de la irregularidad de la estrella, IIS0, cuyo nodo central tiene de coordenadas (x0, y0). Se adopta como indicador de irregularidad de la malla (IIC) el valor mínimo de los indicadores de irregularidad obtenidos para cada una de las estrellas de la malla.
El control sobre la irregularidad de la malla es importante, ya que al disminuir el IIC, disminuye el paso del tiempo. Sin embargo, mediante un proceso adaptativo de la malla (an˜adiendo nodos de forma selectiva [35]) se pueden conseguir índices de irregularidad más próximos a la unidad y, por tanto, pasos de tiempo mayores y menores relaciones de dispersión.
5Resultados numéricos5.1Irregularidad y estabilidadEn esta subsection se muestran los resultados obtenidos al resolver el problema definido en 2.1 en Ω=[0, 1]×[0, 1]⊂R2, con las condiciones de contorno Dirichlet y las condiciones iniciales
utilizando irregulares (ver figuras 4 y 5) de 121 nodos. La solución analíticaLa función de ponderación utilizada ha sidoy el criterio de selección de los nodos es el del cuadrante. El error global ha sido calculado para cada paso de tiempo utilizando la siguiente normadonde sol(j) es el valor de la solución aproximada en el nodo j, exac(j) es el valor de la solución exacta en el nodo j, exacmax es el máximo valor de la solución exacta en los nodos interiores de de la malla considerada y NT es el número de nodos del interior.
Los valores de los errores globales, para n=500, al variar ▵t, utilizando la malla irregular de 121 nodos de la figura 4, con IIC=0, 6524, se muestran en la tabla 1.
Los valores de los errores globales, para n=500, al variar ▵t, utilizando la malla irregular de 121 nodos de la figura 5, con IIC=0, 8944, se muestran en la tabla 2.
5.2Resultados para distintas longitudes de ondaLa tabla 3 muestra los resultados de los errores globales de U y V, para n=500 y ▵t=0, 01, al resolver el problema definido en 2.1, utilizando una malla regular de 121 nodos (figura 3), para funciones cuyas soluciones analíticas son: 39 y
Errores globales para distintas soluciones analíticas del problema 2.1, n=50;▵t=0, 01
| Ecuación (sol. anal.) | Error global U | Error global V |
| 42 | 1, 646×10−6 | 1, 778×10−6 |
| 39 | 1, 232×10−5 | 2, 443×10−5 |
| 43 | 4, 081×10−4 | 2, 001×10−4 |
| 44 | 9, 188×10−2 | 8, 647×10−2 |
| 45 | 3, 035×10−1 | 3, 275×10−1 |
| 46 | 4, 942×10−1 | 4, 999×10−1 |
Los valores obtenidos de la dispersión en cada una de las estrellas de la malla irregular con IIC=0, 8944 de las ondas P, para tres valores distintos de del ángulo φ=0, 45, 90, se han representado sobre los nodos en la figura 6. Puesto que los valores en los tres casos únicamente difieren a partir de la tercera cifra decimal la figura 6 solo muestra el resultado para el caso φ=45
6Conclusiones
En este artículo se muestran las expresiones explíctas, en diferencias finitas generalizadas, para la propogación de ondas sísmicas en 2-D. Se estudia la estabilidad, obteniéndose un criterio en función de la razón de velocidades y los coeficientes de las fórmulas explícitas.
Se analiza la dispersión y se relaciona con la irregularidad de la estrella utilizando el índice de irregularidad de la malla. La utilización de mallas irregulares, según la geometría del problema, puede suponer dispersiones elevadas en algunas estrellas, lo cual va asociado con valores pequen˜os del IIC y un menor paso de tiempo. En este caso se puede aumentar el IIC, y por tanto, el paso del tiempo y disminuir las dispersiones, realizando una redefinición de la malla, para lo que se puede utilizar un proceso adaptativo como el presentado por los autores en [35].
Por último es fácil ver como el error aumenta al disminuir el número de puntos de la discretización por longitud de onda, apreciándose que hasta el orden de diez puntos por longitud de onda se obtienen errores generalmente admisibles inferiores al 10%.
Los autores agradecen la ayuda recibida del Ministerio de Ciencia e Innovación de Espan˜a en el proyecto TISMANCA, Ref.: CGL2008-01757/CLI.