Simulación de ondas sísmicas 2D en un medio isótropo, homogéneo y perfectamente elástico

Trabajo de Fin de Grado Simulaci´on de ondas s´ısmicas 2D en un medio is´otropo, homog´eneo y perfectamente el´astico Autor: Miguel Quevedo Mart´ınez Tutora: Ana M. Portillo de la Fuente 9 de julio de 2021 Resumen Este trabajo trata sobre la aplicaci´on de m´etodos num´ericos para la resoluci´on de un problema de dos ecuaciones en derivadas parciales (EDP) en dos dimensiones espaciales con condici´on inicial y en la frontera de un modelo matem´atico de ondas s´ısmicas. Bajo los supuestos de medio homog´eneo e is´otropo, cuyos materiales se mantienen siempre dentro de sus l´ımites el´asticos, se simulan en dos dimensiones, tanto el problema homog´eneo como con t´ermino fuente. En el primer caso, las ondas s´ısmicas se simulan desde las propias condiciones iniciales y en el segundo, se hace con los t´erminos fuente de tipo Ricker y Ormsby. Adem´as, se ha implementado la posibilidad de simular en la misma regi´on varias capas de diferentes materiales. Se plantean los principios f´ısicos que afectan al medio de propagaci´on de las ondas, a continuaci´on se discretiza el modelo, primero en espacio con diferencias finitas en la direcci´on x, la direcci´on zy la derivada cruzada. Despu´es, el problema semidiscreto resultante se reescribe como un problema de primer orden en tiempo y se resuelve con el m´etodo de Strang. Se implementa en entorno Matlab para finalmente llevar a cabo experimentos de simulaci´on, obteniendo de todo ello resultados y conclusiones. Palabras clave Ondas s´ısmicas, mec´anica de los medios continuos, diferencias finitas, m´etodo de splitting, simulaci´on 2D. Abstract This project involves numerical methods which are used to solve a problem of two partial derivative equations (PDE) in two spatial dimensions with initial and boundary conditions on a mathematical model of seismic waves. Assuming that the medium is homogeneous, isotropic and perfectly elastic, both homogeneous and nonhomogeneous seismic problems of waves are simulated in two dimensions. On the first case, this is made starting from the initial conditions, whereas on the second one this is achieved by using Ricker and Ormsby wavelets as source terms. A way to simulate more than one different material layer into the same study region has also been implemented. First of all, the physical principles which affect the propagation media of seismic waves are raised. Then, the model is discretizated in space with finite differences in the xdirection, in the zdirection and in the mixed derivative. The resulting semidiscrete problem is rewritten as a first order problem in time, and solved using the Strang method. Finally, simulation experiments are carried out in the Matlab environment, leading to results and conclusions. Key words Seismic waves, mechanics of continuous media, finite differences, splitting method, 2D simulation. 1 ´ Indice general 1. Introducci´on 6 1.1. Objetivos del trabajo ........................... 8 1.2. Problema matem´atico y distribuci´on del trabajo ........... 8 1.3. Tipos de ondas s´ısmicas ......................... 9 1.3.1. Ondas internas .......................... 9 1.3.2. Ondas superficiales ........................ 11 1.4. Caracter´ısticas del medio ........................ 11 1.4.1. Ley de Hooke ........................... 13 1.4.2. Par´ametros de Lam´e ....................... 15 1.4.3. Ejemplos de materiales ..................... 17 1.5. Velocidad de propagaci´on ........................ 18 2. Modelo general 19 2.1. Discretizaci´on espacial .......................... 21 2.2. Discretizaci´on temporal ......................... 24 2.3. Estabilidad del modelo .......................... 27 3. Modelo homog´eneo 28 3.1. Condici´on inicial ............................. 28 3.2. Experimentos ............................... 28 3.2.1. Experimentos iniciales ...................... 28 3.2.2. Simulaci´on de materiales reales ................. 34 4. Modelo no homog´eneo 49 4.1. Integraci´on num´erica: Regla de Gauss de tres nodos ......... 49 4.2. Tipos de t´ermino fuente ......................... 49 4.2.1. Ricker ............................... 49 4.2.2. Ormsby .............................. 52 5. Conclusiones 54 5.1. Repercusiones ............................... 54 5.2. L´ınea futura ................................ 54 Bibliograf´ıa 56 2 ´ Indice de figuras 1.1. Representaci´on del terremoto y posterior maremoto de Lisboa en 1755. Autor desconocido. ........................ 7 1.2. Distintos tipos de ondas y la forma de propagaci´on de las mismas a trav´es de las capas interiores de la Tierra. [1]............. 10 1.3. Representaci´on de la propagaci´on de una onda primaria. [5]..... 10 1.4. Representaci´on de la propagaci´on de una onda secundaria. [5]. . . . 11 1.5. Representaci´on de la propagaci´on de una onda Love (a) y una onda Rayleigh (b). ............................... 11 1.6. Representaci´on de la fuerza aplicada en un elemento infinitesimal de un s´olido, la tensi´on resultante y su descomposici´on en componente normal y tangencial. ........................... 12 1.7. Estado tensional de un elemento infinitesimal en el interior de un cuerpo. ................................... 12 1.8. Diagrama del ensayo a tracci´on de un material. ............ 14 2.1. Regi´on espacial o dominio R....................... 21 2.2. Nodo jl y sus adyacentes. ........................ 22 3.1. Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con densidad ρ= 2500 kg/m3,λ= 20 GPa y µ= 20 GPa en el instante t= 0.400 s.................. 29 3.2. Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con ρ= 5000 kg/m3,λ= 20 GPa yµ= 20 GPa en el instante t= 0.400 s......................... 30 3.3. Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con ρ= 2500 kg/m3,λ= 100 GPa yµ= 20 GPa en el instante t= 0.400 s......................... 31 3.4. Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con ρ= 2500 kg/m3,λ= 20 GPa yµ= 100 GPa en el instante t= 0.400 s......................... 32 3.5. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de varias capas en el instante t= 0.400 s............ 33 3.6. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de varias capas en el instante t= 0.400 s............ 34 3.7. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de aceite en el instante t= 0.250 s................ 35 3.8. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de agua en el instante t= 0.250 s................. 36 3.9. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de aceite en el instante t= 1.5s................. 36 3 3.10. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de agua en el instante t= 1.5s.................. 37 3.11. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de anhidrita en el instante t= 0.250 s.............. 38 3.12. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de arenisca en el instante t= 0.250 s............... 39 3.13. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de calcita en el instante t= 0.250 s................ 40 3.14. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de caliza en el instante t= 0.250 s................ 41 3.15. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de cuarzo en el instante t= 0.250 s................ 42 3.16. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de dolomita en el instante t= 0.250 s.............. 43 3.17. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de esquisto en el instante t= 0.250 s............... 44 3.18. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de feldespato en el instante t= 0.250 s.............. 45 3.19. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de pirita en el instante t= 0.250 s................ 46 3.20. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de plagioclasa en el instante t= 0.250 s............. 47 3.21. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de siderita en el instante t= 0.250 s............... 48 4.1. Forma de la ond´ıcula de tipo Ricker (amplitud frente al tiempo) y filtro de frecuencias para un f0= 25 Hz. [17]............. 50 4.2. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ricker de f0= 20 Hz en t= 0.300 s.............. 51 4.3. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ricker de f0= 25 Hz en t= 0.300 s.............. 51 4.4. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ricker de f0= 35 Hz en t= 0.300 s.............. 52 4.5. Forma de la ond´ıcula de tipo Ormsby (amplitud frente al tiempo) y filtro de frecuencias para la combinaci´on de frecuencias 5-10-40-45 Hz, [17]..................................... 53 4.6. Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ormsby de f= 5,10,40,45 Hz en t= 0.300 s......... 53 4 ´ Indice de tablas 1.1. Propiedades f´ısicas de materiales frecuentes en el interior terrestre . 17 1.2. Velocidades a trav´es de distintos materiales .............. 18 2.1. Valores de αij para cada material de ejemplo. ............. 20 2.2. Relaci´on de tama˜nos de paso k hque aseguran la estabilidad por material. 27 3.1. Distribuci´on de capas para la primera prueba. ............. 29 3.2. Distribuci´on de capas para la segunda prueba. ............. 30 3.3. Distribuci´on de capas para la tercera prueba. ............. 30 3.4. Distribuci´on de capas para la cuarta prueba. ............. 31 3.5. Distribuci´on de capas para la quinta prueba. ............. 33 3.6. Distribuci´on de capas para la sexta prueba. .............. 34 3.7. Distribuci´on de capas para la prueba con aceite. ........... 35 3.8. Distribuci´on de capas para la prueba con agua. ............ 35 3.9. Distribuci´on de capas para la prueba con anhidrita. .......... 37 3.10. Distribuci´on de capas para la prueba con arenisca. .......... 38 3.11. Distribuci´on de capas para la prueba con calcita. ........... 39 3.12. Distribuci´on de capas para la prueba con caliza. ............ 40 3.13. Distribuci´on de capas para la prueba con cuarzo. ........... 41 3.14. Distribuci´on de capas para la prueba con dolomita. .......... 42 3.15. Distribuci´on de capas para la prueba con esquisto. .......... 43 3.16. Distribuci´on de capas para la prueba con feldespato. ......... 44 3.17. Distribuci´on de capas para la prueba con pirita. ............ 45 3.18. Distribuci´on de capas para la prueba con plagioclasa. ......... 46 3.19. Distribuci´on de capas para la prueba con siderita. ........... 47 4.1. Distribuci´on de capas para todos los experimentos con la fuente de tipo Ricker. ................................ 51 5 1. Introducci´on El t´ermino “onda s´ısmica1” suele asociarse con terremotos o se´ısmos, que las generan de forma natural. No obstante, adem´as de las naturales, existen causas artificiales que pueden provocarlas. En el caso de los terremotos, las ondas s´ısmicas son generadas por los movimientos de las placas litosf´ericas, formadas por la corteza y la parte m´as superficial del manto terrestre. A una velocidad de unos cent´ımetros al a˜no, dichas placas litosf´ericas se mueven sobre una capa fluida denominada astenosfera. Una velocidad aparentemente lenta, pero suficiente para provocar rupturas s´ubitas en las placas, que se hallan en constante movimiento, donde estas se acercan (terremotos profundos2), rozan entre ellas (terremotos cercanos a la superficie) o se alejan entre s´ı (formando las dorsales oce´anicas) [1]. En la pen´ınsula ib´erica, la cercan´ıa a la zona de contacto entre las placas europea y africana lleva a que se produzcan terremotos desde el mar de Albor´an hasta el cabo de San Vicente (Lisboa, 1755, figura 1.1). La cantidad de terremotos que ocurren en Espa˜na alcanza una media de 3000 al a˜no3, aunque la mayor´ıa de estos son de baja intensidad. Las zonas de mayor actividad s´ısmica en la pen´ınsula son: los Pirineos y la costa sur y sureste del Mediterr´aneo (Lorca, 2011). Adem´as, en las Islas Canarias, de origen volc´anico, tambi´en tiene lugar una gran actividad s´ısmica. Existe una gran variedad de fen´onemos que pueden dar lugar a la aparici´on de ondas s´ısmicas y que nada tienen que ver con los movimientos naturales de las placas en la litosfera. Si el origen de estas ondas es artificial, es decir, producto de la actividad humana4, resulta a´un m´as interesante estudiar y poder simular las consecuencias de dicha actividad. Algunos ejemplos que causan esta “sismicidad inducida5” son: [3] Los lagos artificiales. 1Onda s´ısmica: “tipo de onda el´astica fuerte en la propagaci´on de perturbaciones temporales del campo de tensiones que generan peque˜nos movimientos en las placas tect´onicas.” [1] 2Esta clase de movimiento provoca los efectos m´as profundos (en ocasiones hasta 600km de profundidad). Cuando las placas se acercan y chocan frontalmente, se deforman y se superponen la una sobre la otra. Jap´on en 2011, Sumatra en 2004 o Chile en 1960 son ejemplos de los terremotos m´as devastadores, todos provocados por esta clase de movimiento. El caso extremo es en el que dos continentes chocan, formando entonces cadenas monta˜nosas (como la cordillera del Himalaya entre la India y Asia, con terremoto en Nepal en 2015). 3Existe un recurso en la web del Instituto Geogr´afico Nacional mediante el cual se pueden visualizar los terremotos en la pen´ınsula ib´erica y zonas pr´oximas de hasta los ´ultimos 30 d´ıas: ign.es/web/resources/sismologia/tproximos/prox.html 4La s´ısmica es la rama de la sismolog´ıa dedicada al estudio de estas ondas artificiales. 5Traducido del t´ermino ingl´es induced seismicity [3], empleado para referirse a los terremotos menores causados por la actividad humana. 6 La miner´ıa. Los pozos de residuos. La extracci´on y el almacenamiento de hidrocarburos. La extracci´on de agua subterr´anea. La extracci´on de energ´ıa geot´ermica. La fractura hidr´aulica. Las explosiones y los ensayos nucleares. La actividad humana en las urbes con [4]:  El tr´afico rodado.  Los trenes de metro o convencionales.  Algunos eventos multitudinarios (footquake6).  Los conciertos de m´usica.  Los espect´aculos pirot´enicos (fuegos artificiales). Figura 1.1: Representaci´on del terremoto y posterior maremoto de Lisboa en 1755. Autor desconocido. Ante un fen´omeno tan notable y frecuente, queda patente la necesidad e importancia de herramientas fiables y precisas para la predicci´on del comportamiento de los suelos. A lo largo de este trabajo, se abordar´an diferentes experimentos de simulaci´on que ser´an comparados (siempre que sea posible) con materiales reales. Podr´a entonces comprobarse si tanto el modelo como los m´etodos de resoluci´on del mismo se ajustan a la realidad y si, por tanto, pueden hacerse predicciones basadas en ellos. 6Traducido del ingl´es, vendr´ıa a significar “terremoto de pisadas”, que puede ocurrir si se da una concentraci´on elevada de personas y se suceden saltos o pasos acompasados. 7 longitudinal, y es propio de cada material. Figura 1.8: Diagrama del ensayo a tracci´on de un material. A partir del l´ımite de proporcionalidad, al retirar la carga, el material recupera la forma siguiendo la misma pendiente que en el tramo (0, σp). En el diagrama se han representado en azul tres ejemplos de deformaci´on permanente, partiendo desde casos de una carga aplicada (A, B y C) hasta que se retira por completo (A’, B’ y C’), siempre manteniendo la misma pendiente o ´angulo α. En estos casos se dice que el material ha plastificado. Una vez que se alcanza el valor de tensi´on m´axima del material, σm´ax, el material sufre un estrechamiento por la secci´on por la que, seguidamente, se rompe (punto r del diagrama). Una onda s´ısmica o cualquier onda mec´anica que se propaga a trav´es de un material produce esfuerzos tanto de tracci´on como de compresi´on12. Los par´ametros caracter´ısticos de los materiales se obtienen de forma experimental a trav´es de ensayos de tracci´on en los que se supone que el material tiene un comportamiento similar a compresi´on y que dichos par´ametros son los mismos, pero esto no siempre es as´ı. Si en la figura 1.8 se representase la curva debida al efecto de la compresi´on en el material, esta se situar´ıa en la zona negativa de ambos ejes, con la misma pendiente E, para tensiones y deformaciones negativas. Existen materiales, sobre todo de car´acter fr´agil, cuyo l´ımite el´astico es diferente en funci´on de si trabajan a tracci´on o a compresi´on (σt6=σc). Para este trabajo se considera que estos l´ımites son iguales o que, de ser distintos, nunca se alcanza el menor de ellos en valor 12Por convenio, las tensiones debidas a la compresi´on son negativas: σ < 0. Ocurre de igual manera con las deformaciones derivadas de cada tipo de esfuerzo, siendo positivas para tracci´on (el material se estira) y negativas para la compresi´on (el material se encoge). 14 absoluto. 1.4.2. Par´ametros de Lam´e Aplicando una tensi´on normal a un cuerpo el´astico en cada una de las tres direcciones del triedro, dicho cuerpo experimenta una deformaci´on en todos los ejes. [6] Tensi´on en el eje x: εx=σx Eεy=−νσx Eεz=−νσx E Tensi´on en el eje y: εy=σy Eεx=−νσy Eεz=−νσy E Tensi´on en el eje z: εz=σz Eεx=−νσz Eεy=−νσz E 15 En cada caso, la deformaci´on coincidente con la direcci´on de la tensi´on aplicada es lineal y depende ´unicamente del valor de la tensi´on y de E. Para las otras dos direcciones, depende tanto de la tensi´on y de Ecomo de νo coeficiente de Poisson. Esta constante, que se obtiene de forma experimental, expresa cu´anto var´ıa la secci´on normal a la direcci´on de la tensi´on aplicada. Para materiales is´otropos, como los casos que se ver´an, este valor es igual en todas las direcciones del plano normal aσi. Sumando los valores anteriores para el caso general de un estado tensional en un sistema de ejes cartesianos, se obtienen las leyes generalizadas de Hooke: εx=1 E[σx−ν(σy+σz)],(1.6) εy=1 E[σy−ν(σx+σz)],(1.7) εz=1 E[σz−ν(σx+σy)].(1.8) Asimismo, existe una relaci´on entre las tensiones tangenciales y los ´angulos de deformaci´on, que es la siguiente: γxy =τxy µ, γxz =τxz µ, γyz =τyz µ.(1.9) Se introduce el t´ermino µ, que representa el m´odulo de rigidez transversal, tambi´en llamado m´odulo de cizalladura13. En las ecuaciones de Lam´e, se despeja la tensi´on de las ecuaciones generalizadas de Hooke para dejarlas en funci´on de las deformaciones: σx=λe + 2µεx,(1.10) σy=λe + 2µεy,(1.11) σz=λe + 2µεz.(1.12) Como consecuencia, aparecen diferentes t´erminos, como e, coeficiente de variaci´on volum´etrica. El coeficiente ees constante e invariante; es decir, es independiente de la base vectorial utilizada para expresar las tensiones y deformaciones. e=εx+εy+εz(1.13) Los otros dos t´erminos que aparecen en las ecuaciones de Lam´e son los dos par´ametros de Lam´e que dan nombre a este apartado. El primer par´ametro de 13En ocasiones se denota con “G”. 16 Lam´e es λ, y el segundo es el m´odulo de cizalladura ya visto en las ecuaciones (1.9), µ. Ambos son constantes que dependen de las caracter´ısticas mec´anicas de cada material. M´as concretamente, dependen del m´odulo de Young y del coeficiente de Poisson: λ=νE (1 + ν)(1 −2ν),(1.14) µ=E 2(1 + ν).(1.15) 1.4.3. Ejemplos de materiales Para poder caracterizar el comportamiento de las ondas, se necesita conocer algunas de sus propiedades mec´anicas, como los par´ametros de Lam´e. Con el objetivo de hacer esto de la manera m´as realista posible, se ha recurrido a [9], donde se encuentra una recopilaci´on de las propiedades mec´anicas de materiales que aparecen con bastante frecuencia en el subsuelo. Las magnitudes mostradas en la Tabla 1.114, de izquierda a derecha, son: densidad, primer y segundo par´ametros de Lam´e, m´odulo de Young, m´odulo de compresibilidad15 y coeficiente de Poisson. Tabla 1.1: Propiedades f´ısicas de materiales frecuentes en el interior terrestre Material ρ λ µ E K ν [kg/m3] [GP a] [GPa] [GPa] [GPa] [−] Aceite 830 1.60 - - 1.6 0.500 Agua salada 1030 2.39 - - 2.3 0.500 Anhidrita 2980 26.00 29.0 72.0 45.0 0.230 Arenisca 2500 1.00 15.5 68.5 16.5 0.075 Calcita 2710 56.00 32.0 84.0 77.0 0.320 Caliza (10pu) 2540 35.50 17.5 188.5 54.0 0.330 Cuarzo 2650 8.00 44.0 95.0 37.0 0.070 Dolomita 2870 65.00 45.0 117.0 95.0 0.300 Esquisto (5pu) 2500 13.50 10.5 90.0 26.0 0.270 Feldespato 2620 28.00 15.0 40.0 37.5 0.320 Pirita 4930 59.00 132.0 305.0 147.0 0.150 Plagioclasa 2630 59.00 26.0 70.0 76.0 0.350 Siderita 3960 90.00 51.0 135.0 124.0 0.320 En los casos del aceite y del agua salada, ´unicos ejemplos de fluidos en la Tabla 1.1, tanto el m´odulo de cizalladura (2 º par´ametro de Lam´e) y el m´odulo de Young tienen valor nulo. Esto se debe a que los fluidos no soportan esfuerzos cortantes (µ= 0) ni esfuerzos a tracci´on (E= 0). 14Fuente original: [9]. Para aquellos par´ametros en los que en la tabla original se da un intervalo posible de valores, se ha escogido el valor medio. 15Expresa con qu´e facilidad var´ıa el volumen del material al ser comprimido uniformemente. 17 1.5. Velocidad de propagaci´on En [9], tambi´en se incluyen los valores de las velocidades de las ondas P y S por material. Estos valores se calculan con los par´ametros de Lam´e y la densidad: vP=p(λ+ 2µ)/ρ, (1.16) vS=pµ/ρ, (1.17) que corresponden a las velocidades de las ondas P y S, respectivamente. Retomando los materiales de ejemplo, la velocidades de las ondas por cada uno de ellos se muestran en la Tabla 1.2. Tabla 1.2: Velocidades a trav´es de distintos materiales Material vPvS [km/s] [km/s] Aceite 1.226 0 Agua salada 1.507 0 Anhidrita 5.299 3.120 Arenisca 2.500-4.500 1.725-3.103 Calcita 6.645 3.436 Caliza (10pu) 3.800-6.500 1.900-3.250 Cuarzo 6.008 4.075 Dolomita 7.349 3.960 Esquisto (5pu) 1.800-5.000 1.000-2.777 Feldespato 4.685 2.393 Pirita 8.094 5.174 Plagioclasa 6.487 3.144 Siderita 6.963 3.589 Nuevamente se aprecian valores nulos en el caso de la velocidad de las ondas S en los ejemplos de fluidos. Es evidente, por (1.17), que si su 2 º par´ametro de Lam´e es nulo, tambi´en lo ser´a la velocidad de este tipo de ondas. Esto ya se advirti´o en la subsecci´on 1.3.1: estas ondas se disipan en medios fluidos, al no soportar esfuerzos transversales o cortantes. 18 2. Modelo general Dados los supuestos de medio homog´eneo, is´otropo y perfectamente el´astico estudiados en el cap´ıtulo 1, las ecuaciones en dos dimensiones que los satisfacen se pueden expresar de la siguiente forma: ∂2ux ∂t2=α11 ∂2ux ∂x2+α22 ∂2ux ∂z2+ 2α12 ∂2uz ∂x∂z +fx, ∂2uz ∂t2=α11 ∂2uz ∂z2+α22 ∂2uz ∂x2+ 2α12 ∂2ux ∂x∂z +fz. (2.1) Se tienen dos ecuaciones lineales1; una para cada cada componente del plano que se estudia, x(componente horizontal) y z(componente vertical). Las variables uxy uzrepresentan el desplazamiento de la onda horizontal y vertical, respectivamente. En cada ecuaci´on se iguala la derivada parcial de segundo orden en tiempo ∂2 ∂t2de uxyuzcon las derivadas parciales de segundo orden en espacio ∂2 ∂x2y∂2 ∂z2de uxy uz, respectivamente, en sus dos primeros t´erminos a la derecha de la igualdad. El tercer t´ermino, con una derivada cruzada ∂2 ∂x∂z de la ude la otra ecuaci´on, relaciona uxyuzhaciendo que ambas ecuaciones se encuentren acopladas. Es decir: que en la ecuaci´on de la componente xaparezca uz, y viceversa. Los par´ametros αij son constantes propias del medio material por el que se desplazan las ondas. Dependen de los par´ametros de Lam´e y de la densidad (estudiados en la secci´on 1.4.2) y se calculan de la siguiente manera: α11 =λ+ 2µ ρ,(2.2) α22 =µ ρ,(2.3) α12 =λ+µ 2ρ.(2.4) Los αij propios de los materiales de ejemplo de la tabla 1.1 se resumen en la tabla 2.1. Los t´erminos fxyfzse denominan “t´erminos fuente” y se emplean para introducir perturbaciones del sistema, en forma de focos de emisi´on de ondas, como hipocentros2de terremotos. Pueden ser dependientes de cualquiera de las variables 1Los m´etodos aplicados para la resoluci´on de este problema son v´alidos tanto para problemas lineales como no lineales. 2El hipocentro es el lugar en el interior de la corteza terrestre donde originan las ondas s´ısmicas. No debe confundirse con el t´ermino epicentro, que es el primer punto en la superficie donde se detectan estas ondas. 19 Tabla 2.1: Valores de αij para cada material de ejemplo. α11 α22 α12 Material hkm2 s2i hkm2 s2i hkm2 s2i Aceite 1.9277 0 0.9639 Agua salada 2.2330 0 1.1165 Anhidrita 28.1879 9.7315 9.2282 Arenisca 12.8000 6.2000 3.3000 Calcita 44.2804 11.8081 16.2362 Caliza (10pu) 27.7559 6.8898 10.4331 Cuarzo 36.2264 16.6038 9.8113 Dolomita 54.0070 15.6794 19.1638 Esquisto (5pu) 13.8000 4.2000 4.8000 Feldespato 22.1374 5.7252 8.2061 Pirita 65.5172 26.7748 19.3712 Plagioclasa 42.2053 9.8859 16.1597 Siderita 48.4848 12.8788 17.8030 espaciales o del tiempo, {fx, fz}=f(x, z, t). En el caso de un modelo homog´eneo (cap´ıtulo 3), estos t´erminos se anulan y las ondas se simulan a partir de las condiciones iniciales (2.9-2.10). Los casos con t´erminos fuente se estudiar´an en el cap´ıtulo 4. Condiciones frontera. Se entiende por frontera el conjunto de elementos que se encuenra en el borde de la regi´on sobre la que se realiza el estudio. Las condiciones frontera elegidas para todos los ejemplos de este trabajo son peri´odicas. Para n= x, z: un(a, z, t) = un(b, z, t), z ∈[c, d] (2.5) ∂xun(a, z, t) = ∂xun(b, z, t), z ∈[c, d] (2.6) un(x, c, t) = un(x, d, t), x ∈[a, b] (2.7) ∂zun(x, c, t) = ∂zun(x, d, t), x ∈[a, b] (2.8) Esto significa que las propiedades de los nodos de la frontera izquierda coinciden con los de la derecha, y los superiores, con los inferiores. Es decir, cuando una onda alcanza un extremo de la regi´on, se refleja en el extremo opuesto. Aplicado a la realidad, simular hasta tal l´ımite carece de sentido; m´as adelante, se adecuar´an las condiciones de cada experimento para no llegar a ello. Condiciones iniciales. Los valores que se dan a los nodos para el instante inicial dependen del estado que se quiera simular. La situaci´on m´as simple es en la que no existe en toda la regi´on ninguna perturbaci´on que se propague. Para el problema homog´eneo, las perturbaciones o impulsos iniciales que originen las ondas deben introducirse a partir de las condiciones iniciales. En cualquier caso, dichas condiciones iniciales deben satisfacer las condiciones frontera. ux(x, z, 0) = u0(x, z), ∂tux(x, z, 0) = v0(x, z) (2.9) uz(x, z, 0) = w0(x, z), ∂tuz(x, z, 0) = z0(x, z) (2.10) 20 Se trata de un problema de condiciones iniciales y en la frontera. Para resolver num´ericamente estas ecuaciones es necesario librarse de los t´erminos con derivadas parciales. Esto se lleva a cabo discretizando el sistema. Existen diversas formas de hacerlo, pero se ha escogido como modelo el de [10], que discretiza en espacio mediante diferencias finitas y luego en tiempo mediante un m´etodo de splitting. Las condiciones frontera restringen la discretizaci´on espacial, mientras que las condiciones iniciales restringen la discretizaci´on temporal. 2.1. Discretizaci´on espacial Figura 2.1: Regi´on espacial o dominio R. Para estudiar las ondas en una regi´on acotada, se define una regi´on rectangular R, de tama˜no [a, b]×[c, d] (figura 2.1). Dado que Res una regi´on espacial continua con infinitos puntos en su interior, se simplifica su estudio creando una malla de (N+1)×(M+1) puntos o nodos, espaciados entre s´ı con un paso higual para ambos ejes3. Por tanto, fijando N(n´umero de nodos en el eje horizontal) se definen: h=b−a N, M =d−c h En consecuencia, los nodos de la malla se sit´uan en xj=a+ (j−1)h, para j= 1, ..., N + 1 y en zl=c+ (l−1)h, para l= 1, .., M + 1, formando para la variables uxy uzuna matriz de inc´ognitas ujl(t) = u(xj, zl, t). Tambi´en se aplica para cada uno de los t´erminos fuente, que quedan discretizados para los puntos de la malla; fx(xj, zl, t) yfz(xj, zl, t). La aproximaci´on por diferencias finitas requiere una serie de puntos alrededor de ujl, por lo que se forman dos vectores de un,h para n=x, z y paso h. Para cada j= 1, ..., N + 1 se consideran los vectores, de dimensi´on M+ 1: ux,j = (ux(xj, z1), ..., ux(xj, zM+1))T uz,j = (uz(xj, z1), ..., uz(xj, zM+1))T Estirando la matriz de inc´ognitas de dimensi´on (N+1) ×(M+ 1) por columnas, se obtienen vectores de dimensi´on (N+ 1)(M+ 1): ux,h =uT x,1, ..., uT x,N+1T uz,h =uT z,1, ..., uT z,N+1T y haciendo uso de esta misma t´ecnica con el t´ermino fuente, se obtiene el sistema de ecuaciones diferenciales ordinarias de segundo orden en tiempo d2 dt2ux,h uz,h=Aux,h uz,h+fx,h fz,h(2.11) 3Podr´ıan tomarse pasos hxyhzdiferentes entre s´ı, pero se hace de esta forma para simplificar el proceso: hx=hz=h 21 siendo Auna matriz cuadrada de dimensi´on 2(N+ 1)(M+ 1). A su vez, la matriz Ase subdivide en las matrices A1,A2yA3, tambi´en cuadradas y de tama˜no (N+ 1)(M+ 1): A=1 h2A1A2 A2A3(2.12) En [10], se discretiza en espacio mediante dos m´etodos de diferencias finitas: de segundo y de cuarto orden. Ambos procesos llevan a la soluci´on, aunque con distintos grados de error que aqu´ı no se estudiar´an. Estos m´etodos consisten en calcular, para cada nodo de la malla, su valor teniendo en cuenta el de los nodos circundantes. En este trabajo s´olamente se ver´a el m´etodo de diferencias finitas de segundo orden para la discretizaci´on espacial, ya que es totalmente v´alido y m´as sencillo. En el caso de esta discretizaci´on (central de segundo orden), s´olo se utilizan los nodos siguiente y anterior, de la siguiente forma: ∂2ujl ∂x2≈uj−1,l −2ujl +uj+1,l h2 ∂2ujl ∂z2≈uj,l−1−2ujl +uj,l+1 h2 ∂2ujl ∂x∂z ≈uj−1,l−1−uj+1,l−1+uj+1,l+1 −uj −1, l + 1 4h2 o en forma matricial 1 h2  000 1−2 1 000 ,1 h2  010 0−2 0 010 ,1 4h2  −1 0 1 000 1 0 −1 , respectivamente. Figura 2.2: Nodo jl y sus adyacentes. 22 Reescribiendo el sistema de (2.11) con las matrices de (2.12), te tiene: d2ux,h dt2=A1ux,h +A2uz,h +fx,h, d2uz,h dt2=A2ux,h +A3uz,h +fz,h. (2.13) donde puede verse con claridad que s´olamente aparecen ya derivadas temporales. Las matrices A1,A2yA3son matrices c´ıclicas debido a las condiciones frontera. Esto consiste en que en cada fila se repite el mismo bloque de elementos desplaz´andose hacia la derecha a medida que se avanza una fila o, lo que es lo mismo, propag´andose dicho bloque centrado sobre la diagonal de la propia matriz. Concretamente, esta es la manera en la que se construyen las matrices A1,A2yA3: A1=              B1B20··· 0B2 B2B1B20··· 0 0B2B1B20. . . . . ................. . . . . . 0 B2B1B20 0··· 0B2B1B2 B20··· 0B2B1              A2=              0C20··· 0CT 2 CT 20C20··· 0 0CT 20C20. . . . . ................. . . . . . 0 CT 20C20 0··· 0CT 20C2 C20··· 0CT 20              A3=              ˜ B1˜ B20··· 0˜ B2 ˜ B2˜ B1˜ B20··· 0 0˜ B2˜ B1˜ B20. . . . . ................. . . . . . 0 ˜ B2˜ B1˜ B20 0··· 0˜ B2˜ B1˜ B2 ˜ B20··· 0˜ B2˜ B1              donde los valores de B1,B2,˜ B1,˜ B2son: B1=α22C1−2α11IM+1, B2=α11IM+1, ˜ B1=α11C1−2α22IM+1,˜ B2=α22IM+1, 23 Tabla 3.2: Distribuci´on de capas para la segunda prueba. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 - 10 5000 20 20 3.4641 2.0000 Figura 3.2: Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con ρ= 5000 kg/m3,λ= 20 GPa yµ= 20 GP a en el instante t= 0.400 s. A trav´es de un material cuya densidad es el doble y cuyos par´ametros de Lam´e se mantienen respecto al experimento anterior, la onda P y la onda S avanzan m´as lentamente (figura 3.2). En el mismo instante de tiempo que en la primera prueba, el radio que alcanzan las ondas es inferior que en dicha prueba por esta raz´on. De hecho, si se calcula la relaci´on de velocidades entre ambos tipos de onda a trav´es del mismo material, se demuestra que no depende de la densidad: vP vS =qλ+2µ ρ qµ ρ =sλ+ 2µ µ. Es decir, las ondas se retrasan o adelantan seg´un se var´ıa el valor de la densidad, pero las velocidades se mantienen proporcionales entre s´ı. 3.2.1.2. Cambios en los par´ametros de Lam´e Tercera prueba. En esta prueba se cambia el valor de λpor uno cinco veces mayor. Las caracter´ısticas de la capa se muestran en la tabla 3.3. Tabla 3.3: Distribuci´on de capas para la tercera prueba. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 - 10 2500 100 20 7.4833 2.8284 30 Figura 3.3: Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con ρ= 2500 kg/m3,λ= 100 GPa yµ= 20 GP a en el instante t= 0.400 s. A la vista de la figura 3.3, es evidente que distancia entre la onda P y la onda S crece cuando el primer par´ametro de Lam´e (λ) se incrementa. De hecho, la velocidad de la onda S se mantiene igual a la de primera prueba, porque su velocidad no depende de este par´ametro. Puede afirmarse que el primer par´ametro de Lam´e condiciona ´unicamente la velocidad de propagaci´on de las ondas primarias, aumentando la diferencia entre las velocidades de las ondas P y S cuando crece. Cuarta prueba. En este experimento se cambia el valor de µpor uno cinco veces mayor. Las caracter´ısticas de la capa aparecen reflejadas en la tabla 3.4. Tabla 3.4: Distribuci´on de capas para la cuarta prueba. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 - 10 2500 20 100 9.3808 6.3246 31 Figura 3.4: Componentes horizontal y vertical de los vectores desplazamiento a trav´es un material con ρ= 2500 kg/m3,λ= 20 GPa yµ= 100 GP a en el instante t= 0.400 s. El segundo par´ametro de Lam´e (µ) tiene influencia tanto en la velocidad de propagaci´on de las ondas primarias como en la de las secundarias. Sin embargo, esta influencia es diferente en ondas P y S (1.16-1.17). Cuando todos los par´ametros se mantienen y µaumenta (figura 3.4), el radio que alcanzan ondas primarias y secundarias es superior, pero si se comparan las relaciones de velocidad entre ondas P y S de la primera prueba y la actual: primera prueba: rP= 1.9596 km;rS= 1.1314 km;vS vP =2.8284 4.8990 = 57.73 %, cuarta prueba: rP= 3.7523 km;rS= 2.5298 km;vS vP =6.3246 9.3808 = 67.42 %, donde rPyrSson los radios de las ondas primaria y secundaria, respectivamente. Se comprueba que la diferencia relativa entre las velocidades de ambos tipos de onda se reduce al aumentar el segundo par´ametro de Lam´e. 3.2.1.3. Simulaci´on con varias capas Quinta prueba. Es posible simular m´as de una capa si se modifican los par´ametros αij. Durante la discretizaci´on en espacio, al construir las matrices del m´etodo de diferencias finitas A1,A2yA3se introducen los nuevos par´ametros cuando se alcanza la profundidad correspondiente a una nueva capa. Hay que tener en cuenta que al introducir m´as de una capa, sin importar cu´antas, la primera y la ´ultima deben ser del mismo material para satisfacer las condiciones frontera peri´odicas (2.5-2.8). Se utiliza la misma regi´on Rde las anteriores pruebas (10 ×10 km) y se simula la misma cantidad de tiempo (0.4 segundos). Las capas seleccionadas, junto a su profundidad, se muestran en la tabla 3.5. 32 Tabla 3.5: Distribuci´on de capas para la quinta prueba. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 - 4 2500 20 20 4.8990 2.8284 2 - 8 2500 20 100 9.3808 6.3246 3 - 10 2500 20 20 4.8990 2.8284 Figura 3.5: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de varias capas en el instante t= 0.400 s. Se aprecia en la figura 3.5 c´omo la onda surge en el interior de la capa central y se propaga de la misma manera que en la cuarta prueba, pues el material es el mismo. Las otras dos capas son del material de la primera prueba, cuya velocidad de propagaci´on es menore. Como la capa m´as cercana es la superior, la onda primaria la alcanza antes, desencadenando varios fen´omenos. Por un lado, la velocidad se reduce y, aparentemente, el radio aumenta. Sin embargo, esto ´ultimo se debe a que al ser superior la velocidad vPen la capa central, la onda se propaga a lo largo del l´ımite de la capa m´as r´apido de lo que lo hace al cruzarlo y da esa apariencia de mayor radio o alcance. Por otro lado, aparece una onda “reflejo” que rebota en la frontera entre las dos capas. Esta onda rebotada se encontrar´a y se mezclar´a con la onda secundaria y con todas las ondas originales y reflejo que surjan. Otro fen´omeno que se da es el cambio en la intensidad del color entre una capa y otra. El color var´ıa de intensidad seg´un los valores que toman los elementos de las componentes horizontal y vertical del desplazamiento uxyux. Sexta prueba. Ahora se invierte el orden de las capas respecto al del experimento anterior. El resto de condiciones es igual al de la quinta prueba. 33 Tabla 3.6: Distribuci´on de capas para la sexta prueba. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 - 4 2500 20 100 9.3808 6.3246 2 - 8 2500 20 20 4.8990 2.8284 3 - 10 2500 20 100 9.3808 6.3246 Figura 3.6: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de varias capas en el instante t= 0.400 s. A la visa de la figura 3.6, la onda comienza a propagarse por la capa central, cuyo segundo par´ametro de Lam´e es menor que en las capas superior e inferior. Debido a que la velocidad de propagaci´on es m´as lenta en esta capa, la onda tarda m´as en salir de esta. Cuando cruza la frontera, la onda primaria acelera y pierde intensidad en la siguiente capa, mientras se aprecian los mismos reflejos en el lado de la capa central. Cuando cruza la onda secundaria ocurre lo mismo, s´olo que casi ocurre al final de la simulaci´on. En conjunto, habiendo simulado durante el mismo intervalo de tiempo, la propagaci´on final es menor. De hecho, ninguna onda ha alcanzado la capa inferior. 3.2.2. Simulaci´on de materiales reales Los experimentos que se desarrollan en esta subsecci´on, lo hacen en una regi´on R de tama˜no 6 ×6km, dividida en N= 400 nodos con espaciado h= 15 men espacio y con k= 1 ms en tiempo. Se simula durante un intervalo de 0.25 segundos con los materiales de la tabla 1.1. Los minerales son descritos brevemente para explicar sus caracter´ısticas. Aceite. Su car´acter fluido provoca que no soporte esfuerzos cortantes. El segundo par´ametro de Lam´e, µ, es nulo, as´ı como la velocidad de las ondas secundarias. 34 Tabla 3.7: Distribuci´on de capas para la prueba con aceite. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Aceite 6 830 1.60 - 1.226 0 Figura 3.7: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de aceite en el instante t= 0.250 s. El alcance de las ondas en la prueba con aceite es peque˜no (3.7). Se aprecia lo que parecen ondas secundarias, que por definici´on no deber´ıan existir en un fluido. En el siguiente experimento se probar´a otro fluido para ver si se repite este resultado. Agua. El agua de [9] es salada. Al igual que en la prueba con aceite, tanto la velocidad de las ondas secundarias como el valor del m´odulo de cizalladura son nulos. Tabla 3.8: Distribuci´on de capas para la prueba con agua. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Agua sal. 6 1030 2.39 - 2.23 0 35 Figura 3.8: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de agua en el instante t= 0.250 s. En ambos experimentos con fluidos, y a pesar de que las densidades y el primer par´ametro de Lam´e son diferentes, el resultado es muy similar para el intervalo de tiempo que se ha simulado (figuras 3.7 y3.8). Esto, junto a la presencia de lo que parecen ondas secundarias hace pensar que el modelo no es apto para fluidos. Como la velocidad de propagaci´on de las ondas es lenta en el caso de ambos fluidos, se prueba a simular con estos dos materiales un tiempo mayor: hasta el segundo y medio. Figura 3.9: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de aceite en el instante t= 1.5s. 36 Figura 3.10: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de agua en el instante t= 1.5s. La perturbaci´on del hipocentro no se propaga ni pierde intensidad en ninguno de los dos casos (figuras 3.9 y3.10). Esta perturbaci´on puede ser la onda secundaria que no se propaga, pero en cualquier caso este comportamiento carece de sentido. Ahora se aprecia a simple vista que los resultados son diferentes entre los casos del aceite y del agua en cuanto al radio que alcanza la onda P en cada uno. Esto tiene sentido, pero si se calculan los radios para el tiempo de simulaci´on t= 1.5s, se alcanzan en aceite: rP= 1.839 km, en agua: rP= 3.345 km, y no es posible que sea as´ı, ya que en, el caso del aceite, el radio que se muestra en la figura 3.9 es mayor de lo que deber´ıa y, en el caso del agua, la onda de la figura 3.10 permanece en el interior de R, pero seg´un el c´alculo del radio deber´ıa haber sobrepasado ya la frontera en dicho instante. En conclusi´on, el modelo no parece funcionar correctamente con fluidos. El resto de los materiales de la tabla 1.1 es s´olido, por lo que no se espera este problema se reproduzca durante la simulaci´on y que se parezca m´as a los experimentos iniciales con materiales te´oricos. Anhidrita. Esta roca est´a compuesta por el mineral hom´onimo, de f´ormula CaSO4. Pertenece a la clase de los sulfatos. Es incolora o de color blanco, gris o azulado. La anhidrita es com´un en dep´ositos de sal, aunque es extra˜no encontrarla bien cristalizada. Tabla 3.9: Distribuci´on de capas para la prueba con anhidrita. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Anhidrita 6 2980 23 29 5.299 3.120 37 Figura 3.11: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de anhidrita en el instante t= 0.250 s. La forma de las ondas en la figura 3.11 vuelve a parecerse a la de los experimentos iniciales. La densidad tiene un valor superior a la de la primera prueba y es por eso que el radio que alcanzan las ondas es inferior. Por su parte, los par´ametros de Lam´e tienen valores cercanos entre s´ı. Esto hace que la relaci´on entre velocidades de las ondas primaria y secundaria sea: vS vP =3.120 5.299 = 58.88 %. Los radios que alcanzan las ondas en el instante de la figura 3.11 son: rP= 1.3248 km rS= 0.78 km Arenisca. Roca compuesta al menos en un 50 % por elementos detr´ıticos de di´ametros que van de 0.062 a 2 mm, unidos en una matriz detr´ıtica fina. Concretamente, la arenisca utilizada en este experimento tiene 10 p.u.1 Tabla 3.10: Distribuci´on de capas para la prueba con arenisca. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Arenisca 6 2500 1 15.5 3.578 2.490 1Unidades de porosidad: unidad adimensional que mide la porosidad. Establece la relaci´on de volumen de los poros respecto al volumen total del material. Su valor puede ir de 0 a 100, como si se tratase de un porcentaje de porosidad. 38 Figura 3.12: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de arenisca en el instante t= 0.250 s. El valor de λes m´as de 15 veces inferior al de µ, lo que provoca que la relaci´on de velocidades vS vP =2.490 3.578 = 69.59 % sea m´as cercana a 1 que en los casos vistos hasta ahora y que las ondas P y S se encuentren tan pr´oximas (3.12). El valor de µ, aunque superior a λ, no es muy grande, por lo que el radio de alcance de las ondas es peque˜no, en comparaci´on con el del experimento anterior en el mismo instante: rP= 0.8945 km, rS= 0.6225 km. Calcita. Mineral compuesto por carbonato c´alcico, de f´ormula CaCO3. De entre los minerales de carbonato de calcio, este es el m´as estable. Aproximadamente el 4 % en masa de la corteza terrestre es calcita. Tabla 3.11: Distribuci´on de capas para la prueba con calcita. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Calcita 6 2710 56 32 6.645 3.436 39 Figura 3.19: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de pirita en el instante t= 0.250 s. La pirita posee una densidad que casi duplica la que ven´ıan teniendo los materiales anteriores. Los par´ametros primero y segundo de Lam´e son tambi´en elevados en comparaci´on. Destaca µ, que equivale a m´as del doble que λ. Esto se traduce en una reducci´on de la diferencia entre las velocidades de onda primaria y secundaria y en que su relaci´on sea mayor: vS vP =5.174 8.094 = 63.92 %. Grandes valores en los par´ametros de Lam´e producen mayores velocidades, siendo las de la pirita las m´as elevadas de todos los materiales probados. Esto repercute en el valor de radios que alcanzan las ondas en el instante que se muestra en la figura 3.19 son: rP= 2.0235 km, rS= 1.2935 km. Plagioclasa. Grupo de minerales que se diferencian seg´un su porcentaje de albita y de anortita. Pertenece al grupo de los feldespatos. Pueden emplearse en la fabricaci´on de vidrios y de bloques hormig´on. [14] Tabla 3.18: Distribuci´on de capas para la prueba con plagioclasa. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Plagioclasa 6 2630 59 26 6.4966 3.1442 46 Figura 3.20: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de plagioclasa en el instante t= 0.250 s. El valor de la densidad de la plagioclasa es intermedio, mientras que con los valores de λyµocurre lo contrario que con la pirita, ya que el que ten´ıa un valor mayor ahora tiene el menor y viceversa. De hecho, la relaci´on tambi´en se aproxima a 2:1, causando que la diferencia entre las velocidades de onda se reduzca: vS vP =3.1442 6.4966 = 48.4 %, que es la relaci´on m´as baja de entre todos los materiales probados. El alcance de las ondas P y S (3.21) es: rP= 1.6242 km, rS= 0.7861 km. Siderita. Mineral del grupo de la calcita de f´ormula F eCO3. Es el principal mineral de hierro y posee una gran densidad. Tabla 3.19: Distribuci´on de capas para la prueba con siderita. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Siderita 6 3960 90 51 6.963 3.589 47 Figura 3.21: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de siderita en el instante t= 0.250 s. Ha de tenerse en cuenta que la densidad es la mayor de todos los materiales que se han probado en este trabajo. A pesar de ello, el avance de las ondas (figura 3.21) es bastante extenso gracias a los tambi´en altos valores de los par´ametros de Lam´e. Su relaci´on es similar a la de la plagioclasa: rP= 1.7408 km, rS= 0.8973 km. Y la relaci´on entre velocidades de ondas primaria y secundaria: vS vP =3.589 6.963 = 51.54 %. 48 4. Modelo no homog´eneo En este cap´ıtulo las ondas se simulan a trav´es de la fuente, por lo que no hay necesidad de introducirlas desde las condiciones iniciales. En consecuencia, las condiciones en el instante inicial tomar´an valor nulo en todos los puntos. 4.1. Integraci´on num´erica: Regla de Gauss de tres nodos Como ya se anticip´o en la secci´on 2.2, cuando la funci´on primitiva de la fuente no puede expresarse en forma de una suma finita de funciones elementales, es necesario recurrir a un m´etodo de integraci´on num´erica con una regla de cuadratura de un orden de al menos el mismo que el m´etodo de discretizaci´on empleado (splitting de Strang de segundo orden). La regla de Gauss de dos nodos es el m´etodo que se ha elegido para resolver todos estos casos. Esta regla de dos nodos es de orden 4 y su expresi´on general es la siguiente: Zt+k t f(s)ds ≈k 2ft+k 21−1 √3+ft+k 21 + 1 √3 (4.1) Es decir, se eval´ua la funci´on fuente en los puntos indicados en (4.1), donde k es el paso en tiempo, que no var´ıa una vez se elige y cumple con los criterios de estabilidad. Adem´as, tes el instante que se est´e calculando en cada iteraci´on de la simulaci´on. 4.2. Tipos de t´ermino fuente Si se representa el valor de la amplitud de la onda en funci´on del tiempo medida en un punto, como ocurre en los sism´ografos, se observa c´omo esta comienza valiendo cero, despu´es aumenta y finalmente regresa a cero. Matem´aticamente, las funciones que describen estas breves ondulaciones reciben el nombre de “ond´ıculas”. Son similares a una onda. Las expresiones que definen estas ond´ıculas son muy variadas, siendo algunas de las m´as comunes las de Ricker y Ormsby en el campo de la s´ısmica. Los ejemplos citados ser´an los que se estudien en este cap´ıtulo. Una caracter´ıstica com´un a todos ellos es que sus expresiones s´olamente dependen del tiempo. 4.2.1. Ricker De los dos tipos de ond´ıculas que se van a estudiar, las de Ricker son las m´as sencillas por varias razones. Si se observa la forma de la onda en la figura 4.1, existe un pico de amplitud donde se toma referencia del instante en el que “pasa” la onda s´ısmica. Este pico est´a rodeado por otros dos de menor amplitud y de signo 49 contrario al pico central. El resto de ese gr´afico se mantiene a cero. Tiene, por lo tanto, bastante sencilla. Adem´as, la funci´on de esta ond´ıcula (4.2) posee un ´unico par´ametro, f0, que corresponde a la frecuencia del pico de la ond´ıcula. Una onda de Ricker es como sigue: Ricker(t) = (1 −2π2f2 0t2) exp(−π2f2 0t2),(4.2) donde tes el tiempo en segundos. El valor que suele adoptar para la frecuencia f0 es de 25 Hz (ond´ıcula de Ricker de 31 milisegundos), que corresponde al pico del espectro de frecuencias de esta onda, como puede verse en la figura 4.1. Sin embargo, en art´ıculos como [15] o [16] se utilizan las frecuencias de 35 Hz y 20 Hz, respectivamente. En los experimentos con Ricker se utilizar´an las tres frecuencias mencionadas para comparar los resultados que se obtengan. Figura 4.1: Forma de la ond´ıcula de tipo Ricker (amplitud frente al tiempo) y filtro de frecuencias para un f0= 25 Hz. [17] Otro aspecto de su sencillez reside en que su expresi´on (4.2) tiene primitiva. Si se calcula y se aplica la regla de Barrow, se tiene: IRicker(s) = Zt+k t r(s)ds = =Zt+k t1−2π2f2 0s2exp −π2f2 0s2ds =sexp −π2f2 0s2t+k t= = (t+k) exp h−π2f2 0(t+k)2i−texp −π2f2 0t2. (4.3) 4.2.1.1. Diferentes valores para la frecuencia pico Se lleva a cabo la misma simulaci´on, con las capas de la tabla 4.1 varias veces alterando la frecuencia pico en la expresi´on de Ricker. De esta manera, se pretende para poder observar las diferencias al variar el ´unico par´ametro que posee esta expresi´on. 50 Tabla 4.1: Distribuci´on de capas para todos los experimentos con la fuente de tipo Ricker. Capa Material zmax ρ λ µ vPvS [km] [kg/m3] [GP a] [GP a]km/s km/s 1 Esquisto 2.4 2500 14 11 3.7148 2.0494 2 Dolomita 4.8 2870 65 45 7.3489 3.9597 3 Esquisto 6.0 2500 14 11 3.7148 2.0494 Frecuencia de 20 Hz. Frecuencia empleada en [16]. Figura 4.2: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ricker de f0= 20 Hz en t= 0.300 s. Frecuencia de 25 Hz. Frecuencia empleada en [17]. Figura 4.3: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ricker de f0= 25 Hz en t= 0.300 s. 51 Frecuencia de 35 Hz. Frecuencia empleada en [15]. Figura 4.4: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ricker de f0= 35 Hz en t= 0.300 s. La forma de las ondas para el mismo instante y distribuci´on de capas y materiales es id´entica. Los ´unicos cambios que pueden apreciarse aparecen en los valores que toman los elementos de las componentes de los vectores desplazamiento horizontal y vertical. Es en el experimento con f0= 25 Hz cuando dichos elementos alcanzan valores absolutos m´ınimos. Comparando con el problema homog´eneo de 3, los valores que toman los elementos de los vectores desplazamiento son sensiblemente menores en el caso con fuente. Adem´as, las ondas son mucho m´as finas, perturbando menos espacio a su paso por el medio por el que se propagan. 4.2.2. Ormsby Tambi´en posee un pico de amplitud donde se toma el origen de tiempo, rodeado de m´as oscilaciones que Ricker. La expresi´on es totalmente diferente1a (4.2), contando con cuatro par´ametros que la definen en lugar de uno solo (figura 4.5): Ormsby(t) = "(πf4)2 (πf4−πf3)sinc2(πf4t)−(πf3)2 (πf4−πf3)sinc2(πf3t)#+ −"(πf2)2 (πf2−πf1)sinc2(πf2t)−(πf1)2 (πf2−πf1)sinc2(πf1t)#; (4.4) donde f1es la frecuencia de corte bajo, f2es la frecuencia de paso bajo, f3es la frecuencia de paso alto y 1Aparece el t´ermino sinc, funci´on denominada “seno cardinal” y que se define como sinc(x) = sin(x) x. 52 f4es la frecuencia de corte alto. Los valores que suelen adoptar estas frecuencias son de 5, 10, 40 y 45 Hz, respectivamente, como se muestra en la figura 4.5. Figura 4.5: Forma de la ond´ıcula de tipo Ormsby (amplitud frente al tiempo) y filtro de frecuencias para la combinaci´on de frecuencias 5-10-40-45 Hz, [17] La primitiva de esta funci´on no puede expresarse como una suma finita de funciones elementales. Para poder utilizarla en el modelo se le aplica la ya mencionada regla de cuadratura de Gauss de dos nodos. Figura 4.6: Componentes horizontal y vertical de los vectores desplazamiento a trav´es de tres capas de esquisto, dolomita y esquisto con una ond´ıcula de tipo Ormsby de f= 5,10,40,45 Hz en t= 0.300 s. Las capas son las mismas que en los experimentos con la fuente tipo Ricker (tabla 4.1). La principal diferencia con Ricker es que la zona perturbada por el paso de la onda es mayor, comparando el mismo instante. De hecho, es incluso mayor que en el caso homog´eneo considerado. Cabe destacar que el desplazamiento horizontal ux(figura 4.5 izquierda) toma predominantemente valores negativos en los elementos con mayores valores absolutos. En contraste con este hecho se encuentra el desplazamiento vertical uz(figura 4.5 derecha), cuyo espectro es m´as sim´etrico respecto al cero. 53 5. Conclusiones A lo largo de la introducci´on se han planteado los principios que posteriormente se necesita conocer para entender las ecuaciones del modelo y sus restricciones. Se ha descrito el modelo con ecuaciones en derivadas parciales y se ha discretizado para poder hallar una soluci´on num´erica, cosa que se ha llevado a cabo mediante diferencias finitas de segundo orden para la discretizaci´on de las derivadas de segundo orden en espacio y un m´etodo de splitting, tambi´en de segundo orden, para la discretizaci´on en tiempo. Una vez programado el modelo en Matlab, ha sido posible realizar experimentos de simulaci´on tanto para su versi´on como problema homog´eneo como para el no homog´eneo. En ambos casos se ha probado a simular una o m´as capas dentro de la misma regi´on. A partir de los experimentos aplicados al modelo homog´eneo, pueden extraerse ciertas conclusiones relativas al mismo. En primer lugar, el modelo responde de la forma esperada al variar los valores de las magnitudes f´ısicas por las que queda definido. Esto se hizo con materiales te´oricos cuyos valores se asemejaban a los de ciertos materiales reales para, despu´es, probar los resultados ensayando materiales reales. Ha podido comprobarse, al no dar resultados que se correspondan con la realidad, que el modelo no funciona correctamente con materiales fluidos. Sin embargo, es aplicable a materiales s´olidos que existen y son frecuentes en la corteza terrestre. En cuanto al modelo no homog´eneo o con fuente, se ha podido probar otra manera de simular ondas s´ısmicas mediante dos tipos de fuente diferentes. Adem´as, para una de dichas fuentes ha sido necesario aplicar una regla de integraci´on num´erica. 5.1. Repercusiones Entender y tener capacidad de simulaci´on de ondas s´ısmicas, ya sean de origen natural o artificial, es una herramienta ´util para la ingenier´ıa. Comprender el comportamiento de los materiales a su paso para dise˜nar, mantener o reformar cualquier tipo estructura o m´aquina es esencial si se quiere evitar el fallo de estas a causa de los se´ısmos. Es decir: aplicar estos conocimientos en la ingenier´ıa puede evitar p´erdidas de tipo humano, econ´omico, ambiental, etc. 5.2. L´ınea futura Una de las posibles v´ıas de continuaci´on con este trabajo puede ser tanto la introducci´on de nuevos tipos de fuente, como dar mayor complejidad geom´etrica a 54 las capas. Tambi´en se propone dar alg´un tipo de arreglo o soluci´on al modelo para que funcione con materiales fluidos. 55