UNIVERSIDAD CENTROAMERICANA FACULTAD DE CIENCIA TECNOLOGÍA Y AMBIENTE

Análisis del Efecto de No-horizontalidad de Estratos en la Evaluación de Respuesta Sísmica de Sitio, en el Área Noroeste de la Zona Urbana de la Ciudad de Managua, Nicaragua.

Trabajo monográfico para obtener el título de: Ingeniero Civil

AUTORES: García Tapia, Kenneth Javier Hernández Álvarez, Cristopher Amir

TUTOR: PhD. Edwin A. Obando.

Managua, Nicaragua

Febrero, 2014

DEDICATORIAS

Al único y soberano Dios, mi Padre Celestial, mi amigo y mi sustento. Porque me permitió llegar hasta este momento con salud y bienestar, donde confiadamente puedo expresar: Eben-ezer, hasta aquí me ayudó Jehová. Sabiendo que de él emana toda sabiduría, conocimiento e inteligencia la cual fue manifiesta sobre mi vida, para culminar con éxito este gran reto.

A mi madre y amiga, Silvia Tapia. Que con su amor y cuidados me han hecho sentir el hijo más afortunado del mundo. Sabiendo corregir en el momento indicado con toda sutileza y llenando de paz nuestro hogar con su inmensa paciencia y dirección.

A mi papá Trinidad García, que además de ser el hombre que inculcó cuidadosamente cada uno de los valores que hoy rigen mi vida, me enseño a gozar de los frutos merecidos del trabajo propio.

A mis hermanos, porque cada uno ha sabido aportar algo único en cada etapa de mi vida. Siendo fuentes de motivación inmediata en momentos de aflicción.

En general a toda mi familia, porque han creído certeramente que el potencial que Dios puso en mi vida ha sido usado conforme a su voluntad.

Finalmente a todos mis amigos, aquellos que fueron más que motivadores y ayudadores en cada reto de mi vida. Entendiendo claramente lo que Proverbios 17:17 expresa: En todo tiempo ama el amigo, y es como un hermano en tiempo de angustia.

Kenneth Javier García Tapia

I

Al Todopoderoso Elohim por ser la fuente de fortaleza física, mental y espiritual, a Jesucristo por ser la roca sobre la que esta cimentada mi existencia, un fundamento seguro sobre el cual, si los hombres edifican, no caerán y al Espíritu Santo por guiarme e inspirarme.

A mis maravillosas cuatro madres: Marlene, Guissell, Leonor y Marilú Álvarez por ser pilares fundamentales en mi vida, por brindarme todo su apoyo incondicional no solo durante la carrera, sino durante todo el transcurso de mi vida hasta el día de hoy, manifestando compresión, consuelo en momentos de pesar, cuido en la enfermedad y además mostrando su amor mediante la disciplina, cualidad sumamente requerida en este trabajo y en todo lo que habrá de ser el trayecto de esta jornada de la vida.

A mi hermano Nickjail el cual es sin lugar a dudas un hombre admirable un ejemplo de disciplina e intelecto, resultados de su amor a la lectura y dedicación al estudio, mi modelo a seguir.

A mi hermanita Sofía, dotada de una inocencia sin igual, con la cual he compartido toda clase de alegrías y momentos sin pesar, ella ha sido para mí una poderosa fuente motivadora.

Cristopher Amir Hernández Álvarez

II

AGRADECIMIENTOS

Primeramente agradecemos a nuestro creador por habernos permitido a diario el aire que respiramos, el sustento para poder vivir, la sabiduría para decidir y las fuerzas para poder luchar.

También expresamos nuestro profundo agradecimiento al PhD. Edwin Obando, quien de manera desinteresada decidió compartir parte de su conocimiento en el área de Ingeniería Sísmica, teniendo en todo tiempo su disposición y paciencia para abordar cada una de las temáticas referente a este estudio, teniendo la convicción plena que trabajaba con personas anhelantes al conocimiento.

Además, agradecemos de manera general a todos los docentes que dedicaron un tiempo y espacio para brindar el conocimiento teórico y experimental obtenido a lo largo de su carrera profesional, en especial a nuestro Coordinador de la Carrera Ing. Otoniel Baltodano, así como nuestros docentes en: Ciencias básicas como Ing. José María Rodríguez y Elías Rodríguez quienes fundamentaron conocimientos matemáticos básicos. A nuestros docentes especialistas en estructura, carretera, suelo y agua a quienes se les agradece inmensamente.

Kenneth García & Cristopher Hernández

III

TABLA DE CONTENIDO

1.

INTRODUCCIÓN ............................................................................................. 1 1.1. ALCANCES Y LIMITACIONES .................................................................. 3

2.

OBJETIVOS..................................................................................................... 4 2.1. OBJETIVO GENERAL ................................................................................. 4 2.2. OBJETIVOS ESPECÍFICOS ........................................................................ 4

3.

DESCRIPCIÓN DEL ÁREA DE ESTUDIO ...................................................... 5 3.1. DATOS GENERALES DE MANAGUA ....................................................... 5 3.2. TECTÓNICA Y SISMICIDAD ..................................................................... 7 3.3. MARCO GEOLÓGICO-ESTRUCTURAL REGIONAL ................................ 7 3.4. MARCO GEOLÓGICO-ESTRUCTURAL DEL ÁREA DE ESTUDIO .......... 8 3.5. FALLAS GEOLÓGICAS Y LINEAMIENTOS .............................................. 8 3.5.1. Falla de Mateare .................................................................................. 9 3.5.2. Sistema de Nejapa .............................................................................. 9 3.5.3. Sistema de Cofradías .......................................................................... 9 3.5.4. Sistema Meridional .............................................................................. 9 3.6. FALLAS GEOLÓGICAS Y AMENAZA SÍSMICA ........................................ 9 3.7. FALLAS SÍSMICAS Y AMENAZA VOLCÁNICA ...................................... 11 3.8. HISTORIA SÍSMICA Y SISMICIDAD ACTUAL ........................................ 11 3.8.1. Los Terremotos de 1931 y de 1972 ................................................... 11 3.9. SISMICIDAD ACTUAL ............................................................................. 13

4.

REVISIÓN DE LITERATURA ........................................................................ 15 4.1. GENERALIDADES .................................................................................... 15 4.1.1. Causas de los sismos ......................................................................... 15 4.1.2. Tectónica de Placas ............................................................................ 15 4.1.3. Fallas Geológicas ................................................................................ 16 4.1.4. Tipos de fallas ..................................................................................... 17 4.2. ONDAS SÍSMICAS .................................................................................... 17 4.2.1. Ondas de cuerpo ................................................................................. 17 4.2.2. Ondas Superficiales ............................................................................ 18 IV

4.3. RESPUESTA SÍSMICA EN 1D .................................................................. 20 4.3.1. La aproximación lineal equivalente ..................................................... 24 4.4. RESPUESTA SÍSMICA EN 2D .................................................................. 26 4.4.1. Ecuación de Equilibrio ......................................................................... 26 4.4.2. Ecuación Constitutiva .......................................................................... 26 4.4.3. Ecuación de Equilibrio en términos de desplazamiento ...................... 27 4.4.4. Descomposición del campo de desplazamiento ................................. 27 4.4.5. Ecuaciones de ondas desacopladas ................................................... 28 4.4.6. Ondas de cuerpo ................................................................................. 29 4.4.7. Propagación de Ondas en medio anisotrópico .................................... 30 4.4.8. Comparación para casos de una y dos capas .................................... 31 4.4.9. Modelamiento de propagación de ondas ............................................ 32 4.4.9.1. Dominio de Tiempo contra Dominio de Frecuencia ...................... 33 4.4.9.2. Señales reales o sintéticas ........................................................... 33 4.4.9.3. Sismogramas Sintéticos ............................................................... 33 4.4.10. El método de Elemento Finito ........................................................... 34 4.4.10.1. Formulación fuerte ...................................................................... 34 4.4.10.2. Formulación débil ....................................................................... 35 4.4.10.3. Minimización aproximada: Método de Galerkin .......................... 36 4.4.10.4. Discretización del dominio .......................................................... 36 4.4.11. Varios tipos de elementos finitos ....................................................... 36 4.4.12. Algoritmos de integración del tiempo ................................................ 40 4.4.12.1. El método de diferencia central .................................................. 41 4.4.12.2. El método de Newmark .............................................................. 41 4.4.12.3. Método de pasos múltiples ......................................................... 42 4.4.13. Modelado de propagación de ondas en medios infinitos ................... 43 4.4.13.1. Contornos Absorbentes en 2D.................................................... 44 4.4.13.2. Caso bidimensional .................................................................... 45 4.4.14. Efecto de sitio: análisis 1D, 2D .......................................................... 46 4.4.15. Región de absorción ......................................................................... 46 4.4.16. Geometría adaptable y tamaño de la malla ....................................... 47 4.5. ESPECTRO DE RESPUESTA .................................................................. 48 4.5.1. Espectro de Respuesta de Deformación. ............................................ 49 4.5.2. Espectro de Respuesta de Pseudo-Velocidad .................................... 50 V

4.5.3. Espectro de Respuesta de Pseudo- Aceleración ................................ 50 4.5.4. Espectro de Respuesta Combinado D-V-A ......................................... 51 5.

PROCESAMIENTO Y ANÁLISIS DE DATOS ............................................... 53 5.1. DETERMINACIÓN DE MODELOS DE VELOCIDAD ................................ 53 5.2. MODELAMIENTO EN DOS DIMENSIONES ............................................. 55 5.3. CÁLCULO DE RESPUESTA SÍSMICA DE SITIO ..................................... 57 5.3.1. Cálculo de acelerogramas sintéticos en superficie.............................. 58 5.3.2. Proceso de convolución ...................................................................... 58 5.3.3. Cálculo de espectros de respuesta ..................................................... 59

6.

RESULTADOS .............................................................................................. 61 6.1. MODELOS EN 2D ..................................................................................... 61 6.2. FUNCIÓN DE TRANSFERENCIA SÍSMICA .............................................. 66 6.2.1. Validación del uso de elementos finitos para propagación de ondas 1D ....................................................................................................................... 66 6.2.2. Análisis de no-horizontalidad de estratos ............................................ 70 6.3. ACELEROGRAMAS SINTÉTICOS EN SUPERFICIE ............................... 73 6.4. ESPECTROS ELÁSTICOS DE RESPUESTA ........................................... 75

7.

CONCLUSIONES .......................................................................................... 81

8.

RECOMENDACIONES .................................................................................. 83

9.

BIBLIOGRAFÍA ............................................................................................. 84

10.

ANEXOS .................................................................................................... 90

VI

LISTA DE FIGURAS Figura 3.1. Mapeo del área de estudio, mapa de Nicaragua y Managua tomados de INETER. .................................................................................................................. 6 Figura 3.2. Emplazamiento tectónico de Placa Coco y Caribe (W. Strauch, 2000) 7 Figura 3.3. Modelo de terreno del área de Managua y alrededores. Se aprecia la compleja situación geológica estructural de la zona. (W. Strauch, 2000) ............... 8 Figura 3.4. Mapa de zonas de fallas geológicas principales. (Consultants, 1975) 10 Figura 3.5. Zonas de fallas y áreas blancas que requieren futura investigación. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002) ...................... 10 Figura 3.6. Fotografía histórica de la Laguna de Tiscapa. (W. Strauch, 2000) ..... 11 Figura 3.7. Fotos del terremoto de Managua, 1972 colección Steinbrugge. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002) ..................................... 12 Figura 3.8. Acelerograma del terremoto de 1972. (W. Strauch, 2000) ................. 12 Figura 3.9. Epicentros de las réplicas del terremoto de Managua de 1972. (W. Strauch, 2000) ....................................................................................................... 13 Figura 3.10. Red Sísmica Nacional de Nicaragua. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002) ................................................................. 14 Figura 3.11. Vista Tridimensional de la sismicidad de Nicaragua. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002) ..................................... 14 Figura 4.1. Principales zonas tectónicas, lomos oceánicos y zonas de subducción. (Kramer, 1996) ...................................................................................................... 16 Figura 4.2. Falla de San Andrés (falla por desgarramiento). (ROSENBLUETH, 1992) ..................................................................................................................... 16 Figura 4.3. Tipos de fallas según su patrón de ruptura: (a) Fallas normales, (b) Fallas inversas, (c) Fallas transformantes (Kramer, 1996). ................................... 17 Figura 4.4. Deformaciones producidas por las ondas P (Bolt, 1993) .................... 18 Figura 4.5. Deformaciones producidas por las ondas S (Bolt, 1993).................... 18 VII

Figura 4.6. Deformaciones producidas por las ondas L (Bolt, 1993) .................... 19 Figura 4.7. Deformaciones producidas por las ondas R (Bolt, 1993) ................... 19 Figura 4.8. Terremoto de Kermadec de 11 de Junio de 1957 (Kramer, 1996) ..... 19 Figura 4.9. Representación esquemática del modelo tensión-deformación utilizada en el modelo lineal equivalente. (J.F. Semblat & Pecker, 2009). .......................... 20 Figura 4.10. Depósito de suelo estratificado para su estudio mediante el método unidimensional. (Schnabel, 1972) ......................................................................... 21 Figura 4.11. Terminología utilizada en análisis de respuesta de sitio. (J.F. Semblat & Pecker, 2009)..................................................................................................... 24 Figura 4.12. Proceso de iteración para el valor del módulo de corte y la razón de amortiguamiento en el análisis equivalente lineal. (J.F. Semblat & Pecker, 2009) 25 Figura 4.13. Relación entre la velocidad de la onda P y la velocidad de la onda S. .............................................................................................................................. 29 Figura 4.14. Función de transferencia a través de una capa elástica del suelo (caso 1) o a través de dos capas elásticas con diferentes velocidades (casos 2, 3 y 4). (J.F. Semblat & Pecker, 2009) .............................................................................. 31 Figura 4.15. Envolvente del sismograma sintético propuesto por Pousse, et al. (2006) .................................................................................................................... 34 Figura 4.16. Discretización del dominio Ω por elementos finitos Ω𝑒 y una descripción de los desplazamientos nodales (caso 2D). (J.F. Semblat & Pecker, 2009) ......... 37 Figura 4.17. Funciones de forma lineal (arriba) y cuadrática (abajo) elementos finitos. (J.F. Semblat & Pecker, 2009) ................................................................... 38 Figura 4.18. Coordenadas intrínsecas de los elementos finitos cuadriláteros distorsionados: elementos distorsionados (arriba) se definen a partir de elementos de referencia (abajo). (J.F. Semblat & Pecker, 2009). .......................................... 39 Figura 4.19. Varios tipos de elementos finitos de diferente orden p. (Hughes, 1987) .............................................................................................................................. 40 Figura 4.20. Condiciones de estabilidad de integración-tiempo de Newmark (Curnier, 1994; Hughes, 1987). ............................................................................. 43 Figura 4.21. Malla en 2D con amortiguador y resortes en su contorno (Bisch, et al., 1999) ..................................................................................................................... 44

VIII

Figura 4.22. Esquema de elementos infinitos en 2D a lo largo de una dirección (arriba) o de dos direcciones (abajo), propuesto por J.F. Semblat and Pecker (2009). .............................................................................................................................. 45 Figura. 4.23. Amplificación de las ondas sísmicas: fenómenos principales en 1D (izquierda) y estructuras geológicas 2D/3D (derecha). (J.F. Semblat & Pecker, 2009) .............................................................................................................................. 46 Figura 4.24. Esquema del modelo de elementos finitos adaptativo con regiones analizadas y de absorción aplicando el análisis de estructuras de pavimentos. (Ryden & Castaings, 2009) ................................................................................... 48 Figura 4.25. (a) Aceleración del suelo, (b) Respuesta de deformación de tres sistemas simples con 𝜉 = 2% 𝑦 𝑇𝑛 = 0.5; 1; 2 𝑠𝑒𝑔, (c) Espectro de Respuesta de Deformación para 𝜉 = 2% (Chopra, 2007) ............................................................ 49 Figura 4.26. Espectro de Respuesta de Pseudo-Velocidad. (Chopra, 2007) ....... 50 Figura 4.27. Espectro de Respuesta de Pseudo- Aceleración. (Chopra, 2007) ... 51 Figura 4.28. Espectro de Respuesta Combinado D-V-A para el sismo de El Centro, 𝜉 = 2% (Chopra, 2007).......................................................................................... 51 Figura 5.1. Mapa de Ubicación de las zonas de transición seleccionadas ........... 54 Figura 5.2. Modelo de Velocidad de cortante ....................................................... 55 Figura 5.3. Interfaz gráfica de software Comsol Multiphysics para el procesamiento en mecánica de sólidos. ........................................................................................ 56 Figura 5.4. Post-procesamiento en rutinas de Matlab describiendo una función de trasferencia............................................................................................................ 57 Figura 5.5. (a) Rutina en Matlab para cálculo de función de transferencia en condiciones de horizontalidad del suelo (b) función de transferencia generada por algoritmo................................................................................................................ 58 Figura 5.6. (a) Algoritmo en Matlab para generar el acelerograma en superficie (b) función de transferencia sísmica y acelerograma en basamento y en superficie. . 59 Figura 5.7. (Izquierda) Post-procesamiento en Matlab para generar espectros de respuesta (derecha superior) espectro de deformación, (derecha central) espectro de velocidad y (derecha inferior) espectro de aceleración .................................... 60 Figura 6.1. Mapa propuesto por Castillo and Zepeda (2013) mostrando la distribución de zonas. Líneas de transición para creación de perfiles 2D generada por los autores....................................................................................................... 62

IX

Figura 6.2. Perfiles Estratigráficos (2D): (a) Bo. El Bóer, (b) Bo. Jorge Dimitrov, (c) Bo. Recreo Sur. ..................................................................................................... 63 Figura 6.3. Modelo 1, Bo. El Bóer: (superior) Funciones de transferencia sísmica generadas por Comsol y Matlab, (centro) perfil 2D y (inferior) Columnas estratigráficas 1D. ................................................................................................. 67 Figura 6.4. Modelo 3, Bo. Jorge Dimitrov: (superior) Funciones de transferencia generadas por Comsol y Matlab, (centro) perfil 2D y (inferior) Columnas estratigráficas 1D. ................................................................................................. 68 Figura 6.5. Modelo 6, Bo. Recreo Sur: (superior) Funciones de transferencia generadas por Comsol y Matlab, (centro) perfil 2D y (inferior) Columnas estratigráficas 1D. ................................................................................................. 69 Figura 6.6. Perfiles y Funciones de transferencia sísmica más representativos: (superior) Bo. El Bóer, (centro) Bo. Jorge Dimitrov y (inferior) Bo. El Recreo. ...... 71 Figura 6.7. Acelerograma deconvolucionado del terremoto de Managua de 1972 (Castillo & Zepeda, 2013). ..................................................................................... 73 Figura 6.8. Acelerogramas sintéticos representativos durante la transición: (superior) Bo. El Bóer, (centro) Bo. Jorge Dimitrov, (inferior) Bo. Recreo Sur....... 74 Figura 6.9. Espectros Elásticos de Respuesta Bo. El Bóer-Edificio INSS: (a) desplazamiento, (b) velocidad y (c) aceleración.................................................... 78 Figura 6.10. Espectros Elásticos de Respuesta Bo. Jorge Dimitrov: (a) desplazamiento, (b) velocidad y (c) aceleración.................................................... 79 Figura 6.11. Espectros Elásticos de Respuesta Bo. Recreo Sur: (a) desplazamiento, (b) velocidad y (c) aceleración............................................................................... 80

X

LISTA DE TABLAS Tabla 4.1. Relaciones entre los parámetros definidos por el comportamiento elástico lineal isotrópico de los materiales. (J.F. Semblat & Pecker, 2009) ........................ 27 Tabla 4.2. Valores típicos de velocidades para las ondas P y S para diferentes tipos de suelo. (J.F. Semblat & Pecker, 2009) ............................................................... 30 Tabla 6.1. Agrupación de modelos y velocidades promedios para cada estrato .. 65 Tabla 6.2. Resumen de pendientes (%) para cada zona de transición ................. 66 Tabla 6.3. Parámetros geométricos para modelamiento en Comsol .................... 70 Tabla 6.4. Registro de aceleraciones máximas .................................................... 75 Tabla 6.5. Valores máximos de espectros ............................................................ 76

XI

RESUMEN En esta investigación se realiza un análisis de respuesta lineal elástica considerando la presencia de estratos no-horizontales (efecto 2D) que puede presentar el suelo. Inicialmente el análisis se realiza en términos de función de transferencia sísmica, lo cual es comparado con la respuesta sísmica considerando estratos horizontalmente infinito (efecto 1D). El análisis del efecto 2D está basado en el modelamiento de secciones transversales de suelo de la respuesta lineal elástica a través del método de elementos finitos en dominio de frecuencias en un rango de 1 Hz a 10 Hz. Para el área noroeste de la zona urbana de la ciudad de Managua se analizaron 9 secciones transversales que contienen regiones de transición con profundidad de basamento variable. Las funciones de transferencia muestran un comportamiento sumamente distorsionado cuando las pendientes son mayores al 10%, mientras que para valores menores la respuesta es casi 1D. Este efecto también se observó en los acelerogramas y espectros elásticos de respuesta sugiriendo que en caso de transiciones es requerido modelar el efecto 2D, estimando así la respuesta de sitio de forma más realista. ABSTRACT This research is focused on the effect of no-horizontal strata conditions in a ground response analysis (2D effect). Initially the analysis is carried out in terms of seismic transfer function and the resultas are later compared to the seismic ground response assuming no-horizontal condition (1D effect). The 2D analysis is based on cross sections that contains a transition region which is modeled in the frequency domain in a 1 to 10 Hz frequency range. Within the transition region the seismic transfer function depicts a clear distortion when the slope of the transition in higher that 10%, while for lower slope values an almost 1D response occurs. This effect is also evident in the time history accelerograms and in the elastic response spectra which indicates the need of performing 2D ground response modelling when strong lateral variation is present.

XII

1. INTRODUCCIÓN

Nicaragua es un país, que por su ubicación dentro del cinturón de fuego del Pacífico, está expuesto al devastador efecto de los terremotos. La actividad de subducción está dominada por el acomodamiento de las placas coco y caribe, provocando la liberación de energía que se propagan a modo de ondas desde la fuente hasta la superficie causando graves daños en las construcciones. Localmente la ciudad de Managua capital de Nicaragua se ubica sobre un sistema de fallas superficiales con un gran potencial de generar terremotos destructivos. Uno de estos terremotos fue el ocurrido el 23 de Diciembre del año 1972 en la ciudad de Managua teniendo una magnitud de 6.2 Ms en la escala Richter. (W. Strauch, 2000). Después del terremoto de 1972 se han llevado a cabo varios estudios en la ciudad de Managua relacionados al comportamiento dinámico de los suelos durante eventos sísmicos. Faccioli, Santoyo, and León (1973), propusieron caracterizaciones de sitio en base a registros SPT. Además se han hecho propuestas generales de distribución de zonas en el área urbana noroeste de Managua en base a modelos de velocidad de corte (Castillo & Zepeda, 2013); sin embargo todos los estudios realizados hasta la fecha se han fundamentado en el análisis de respuesta sísmica de sitio usando el análisis en una dimensión asumiendo que los estratos son horizontalmente infinitos. En sitios donde la nohorizontalidad de los estratos es marcada, también conocida como efecto 2D, los métodos tradicionales basados en la teoría 1D pueden no ser tan efectivos causando una estimación errónea de la respuesta del suelo. De la misma manera estudios realizados por Obando Hernández (2011) sugieren que los suelos de la Ciudad de Managua no pueden ser considerados homogéneos, debido que la profundidad de la roca elástica o impedancia entre dos estratos se encuentran a profundidades mayores a los 10 m. Como soporte de que puedan existir discrepancias significativas en las estimaciones de las respuestas sísmicas de sitio, se tienen estudios previos elaborados en otros países. Por ejemplo en Venezuela J. J. Hernández, Schmitz, Delavaud, Cadet, and Domínguez (2011), detallan que estudios en 2D proyectan datos más fidedignos de las respuestas de perfiles estratigráficos, además en la zona de la Cerdaña entre España y Francia, el estudio de Tapia, Macau, and Figueras (2007) demostró que existen diferencias entre las simulaciones 1D y 2D debidas principalmente a la geometría del terreno. Dada las conclusiones de dichos estudios aplicando el modelamiento en dos dimensiones, se hace necesario un estudio que demuestre la existencia de posibles variaciones entre las simulaciones en una y dos dimensiones. 1

En lo referente a estudios de comportamiento dinámicos del suelo, es requerido determinar de forma más realista las aceleraciones espectrales, el cual es un dato fundamental para realizar diseños sísmicos de estructuras tomando en cuenta los parámetros de vibración (nivel de amplificación y periodos naturales de vibrar) respecto al basamento rocoso, evitando la sobre y sub-estimación del impacto sísmico. Para este estudio se pretende analizar el efecto de no-horizontalidad de estratos, en el área noroeste de la zona urbana de la capital debido a datos más altos en intensidad para el terremoto de 1972, que datan de la estación ubicada en la refinería ESSO durante la eventualidad sísmica. Para ello se obtendrán 9 secciones transversales (perfiles estratigráficos 2D) basados en los modelos de velocidad de cortante obtenidos por Castillo and Zepeda (2013) en las zonas de análisis definidas por los mismos, para posteriormente definir líneas de transición y observar las posibles variaciones en cada estrato del suelo. Luego a partir de estos perfiles se construirán modelos diferenciados para cada línea de transición, siendo analizados a través del modelamiento en elemento finito en dominio de frecuencia aplicando el Software Comsol Multiphysics y específicamente el módulo de mecánica de sólidos. Los resultados del modelamiento serán presentados en funciones de transferencia sísmica los cuales serán comparados con la respuesta equivalente obtenida a partir del análisis en una dimensión. Seguidamente se utilizará el registro sísmico obtenido en la Refinería ESSO para el terremoto de 1972, como señal de entrada para estimar los espectros elásticos de respuesta tomando en cuenta la variabilidad de estratos y notar si la ocurrencia de la no-horizontalidad en los estratos influye en cambios significativos tanto en las máximas aceleraciones espectrales así como en el ancho del Plateau del espectro.

2

1.1.

ALCANCES Y LIMITACIONES

Con el siguiente trabajo Monográfico se planea generar nueva información referente a respuesta dinámica de los suelos de Managua. Este trabajo está basado en los resultados obtenidos de estudios de microzonificación sísmica previos Castillo and Zepeda (2013) que a su vez se fundamentan en los trabajos de (Escorcia & Ochoa, 2013; Faccioli, et al., 1973; R. O. Hernández, 2009; Obando Hernández, 2011), el cual permitió obtener perfiles 2D a partir de los modelos de velocidades. Estos perfiles se modelan con software fundamentados en procesamiento con elementos finitos (Comsol Multiphysics); además, se establece procedimientos que definen una metodología de trabajo práctica para posteriores estudios en esta área en específico. Así mismo, el estudio tendrá una comparación entre dos modelamientos (1D y 2D) los cuales definirán por acertada o no la hipótesis planteada, determinando con exactitud si en la zona de estudio es meritorio más análisis en dos dimensiones y si es requerido para futuros trabajos expandirla a otras zonas de sensibilidad sísmica. Por otro lado, el modelamiento de los perfiles 2D está basado en la microzonificación y modelos de velocidad de corte promedios propuestos por Castillo and Zepeda (2013), lo que implica que estos datos no fueron obtenidos por visitas de campo aplicada por los autores debido a factores de costo y tiempo. Finalmente es importante mencionar que para el análisis se asume respuesta lineal elástica por lo que no se toma en consideración posible efecto de deformaciones los cual debería ser considerado en futuros estudios.

3

2. OBJETIVOS

2.1. OBJETIVO GENERAL Análisis del Efecto de No-horizontalidad de Estratos en la Evaluación de Respuesta Sísmica de Sitio, a partir de Modelamiento de Propagación de Ondas en Dos Dimensiones, en el Área Noroeste de la Zona Urbana de la Ciudad de Managua, Nicaragua. 2.2. OBJETIVOS ESPECÍFICOS 

Proponer modelos de velocidad de cortante en 2D con variación marca en los estratos basados en estudios previos (Castillo & Zepeda, 2013).



Estimar el grado de aproximación de la respuesta lineal elástica en modelos de 1D para validar la herramienta Comsol sustentada en elementos finitos.



Analizar posibles diferencias en amplificación y contenido de frecuencia usando el análisis en 1 y 2 dimensiones.



Calcular máximas aceleraciones en superficie usando el terremoto de 1972 como registro sísmico de entrada.



Analizar posibles efectos de la no-horizontalidad del suelo en el cálculo de espectros elásticos de respuesta para la zona de estudio.



Proponer un procedimiento para el análisis de amplificación de ondas en dos dimensiones a partir de modelamiento en elemento finito.

4

3. DESCRIPCIÓN DEL ÁREA DE ESTUDIO

3.1.

DATOS GENERALES DE MANAGUA

La ciudad de Managua cabecera del departamento de Managua, fue fundada con el nombre de Villa de Santiago de Managua y es actualmente la capital de la republica de Nicaragua. Su densidad poblacional es de 1, 817,096 habitantes (2004) en su área metropolitana. La ubicación de la Ciudad de Managua está entre las coordenadas 12 o 01´-12o 13´ Latitud Norte y 86o 07´-86o 23´ Latitud Oeste. Sus límites municipales son al Norte con el Lago Xolotlán, al Sur con el Municipio de El Crucero, al Este con el Municipio de Tipitapa, Nindirí y Ticuantepe y al Oeste con el Municipio Villa Carlos Fonseca y Ciudad Sandino. Cuenta con una superficie municipal de 289 km 2 con una superficie urbana de 150.5 km2. La figura 3.1 muestra el mapeo para el área de estudio. Según la Dirección de Planificación/Estadísticas (2007) de la alcaldía de Managua, la capital se encuentra a una altitud mínima de 43 m.s.n.m (metros sobre el nivel de mar) y a una altitud máxima de 700 m.s.n.m. El clima es tropical de sabana, caracterizado por una prolongada estación seca y temperaturas altas durante el año, varían desde 27o C a 34o C.

5

Figura 3.1. Mapeo del área de estudio, mapa de Nicaragua y Managua tomados de INETER.

6

3.2.

TECTÓNICA Y SISMICIDAD

La ocurrencia de terremotos en Nicaragua se debe a que la placa tectónica del Coco choca con la placa tectónica del Caribe, y desciende abruptamente en un ángulo de 80 grados en dirección Noreste bajo el margen pacifico de la placa Caribe. En esta zona de contacto se generan sismos y grandes terremotos con magnitudes hasta 8 Richter. Debajo de Managua, la placa subducida alcanza más de 200 km. En esta profundidad, parte del material de la placa Coco es fundida debido a las altas temperaturas del manto terrestre. (W. Strauch, 2000). El movimiento relativo convergente de estas placas tiene una tasa de 8 centímetros aproximadamente por año. (DeMets, Gordon, Argus, & Stein, 1994), (Figura 3.2).

Figura 3.2. Emplazamiento tectónico de Placa Coco y Caribe (W. Strauch, 2000)

3.3.

MARCO GEOLÓGICO-ESTRUCTURAL REGIONAL

Fisiográficamente, Managua está localizada dentro de la cordillera volcánica y en la porción central de la Depresión o Graben Nicaragüense, un graben poco profundo de más de 300 Km de extensión y 70 km de ancho, con dirección NO-SE, que cruza el territorio nacional en el sector occidental, paralela a la costa del pacifico y a la fosa mesoamericana. Se extiende desde Guatemala hasta el norte de Costa Rica y es rellenado por una espesa secuencia de depósitos volcánicos, volcanoclasticos, aluviales y lacustres (Consultants, 1975). En su parte media se encuentra la cadena volcánica activa de Nicaragua y los grandes lagos: Xolotlán y Cocibolca. 7

3.4.

MARCO GEOLÓGICO-ESTRUCTURAL DEL ÁREA DE ESTUDIO

Managua se ubica dentro de la cordillera volcánica entre los volcanes Apoyeque al noroeste y Masaya al sureste. En ella y sus alrededores se reconocen numerosos montes, volcanes y remanentes de volcanes: Santa Ana, Asososca, Tiscapa, Ticomo, Motastepe, entre otros. El subsuelo de Managua es caracterizado por la presencia de una secuencia volcano-sedimentaria provenientes de los volcanes Masaya, Apoyeque, Apoyo, del lineamiento Miraflores-Nejapa, entre otros. La presencia de numerosos suelos fósiles demuestra la existencia de ciertos periodos de calma entre eventos volcánicos o tectónicos, permitiendo el desarrollo de suelos de varios tipos. (Hradecky, 1997).

Figura 3.3. Modelo de terreno del área de Managua y alrededores. Se aprecia la compleja situación geológica estructural de la zona. (W. Strauch, 2000)

3.5.

FALLAS GEOLÓGICAS Y LINEAMIENTOS

El área de Managua está ubicada dentro de la depresión de Managua (Figura 3.3), esta depresión o graben está limitada por la Falla Cofradía al este y el lineamiento Miraflores-Nejapa al oeste. Hacia el norte el graben se pierde dentro del lago y hacia el sector suroeste el graben es limitado por la Falla Mateare y la Falla Las Nubes, mientras hacia el sur su límite está dentro de las calderas de Las Sierras. La Femina, Dixon, and Strauch (2002) Explican la orientación preferencial NNESSO con desplazamiento lateral izquierdo de las fallas en la cadena volcánica como acomodación de los bloques tectónicos en la cadena volcánica.

8

3.5.1. Falla de Mateare Esta falla de dirección general Noroeste-Sureste, marca el límite occidental de la “Región de Managua” y localmente constituye ese mismo límite en la Depresión de Nicaragua. Es una falla de tipo normal. Las evidencias indican que esta falla tiene su menor desplazamiento en la zona situada al norte de El Crucero (Kuang, 1971), llegando a perderse su traza en esta zona, y que ese desplazamiento aumenta al norte hacia la región de Mateare. 3.5.2. Sistema de Nejapa Una zona de dirección general norte-sur y una anchura promedio de 2,500 metros, se extiende desde el sureste de la Laguna de Apoyeque hasta el sur de la región de Ticomo. Esta zona está caracterizada por un sistema de fracturación el cual da origen a una serie de estructuras volcánicas (Kuang, 1971), que tienen diversas expresiones, siendo más notable aquellas debido al colapso. Este sistema, que ha originado una zona de debilidad estructural por la cual ascendieron soluciones magmáticas a la superficie, desvió la dirección general noroeste-sureste sobre la cual se alinean los volcanes de la Cordillera de los Maribios, en forma paralela a la fosa Mesoamericana. (Arce Velasco, 1973) 3.5.3. Sistema de Cofradías Según Arce Velasco (1973), le ha denominado así a un grupo de fallas dispuestas “en echelon” y que se extiende desde la zona de Tipitapa en dirección sur-suroeste, hasta la región de San Juan de la Concepción. Se trata de fallas normales, cuya dirección de desplazamiento es opuesta al de los anteriores sistemas y que marca de esta manera el límite este de la depresión menor que constituye la “Región de Managua”. 3.5.4. Sistema Meridional Este sistema comprende las fallas normales que limitan por al sur la depresión de la “Región de Managua”, extendiéndose desde la zona al sureste de El Crucero, hasta la región de San Juan y La Concepción. Estas fallas tienen un desplazamiento cuya manifestación topográfica es bastante grande y prácticamente cierran la depresión, con el desplazamiento hacia abajo de los bloques situados hacia el norte de dichas fallas. (Arce Velasco, 1973) 3.6.

FALLAS GEOLÓGICAS Y AMENAZA SÍSMICA

En la ciudad de Managua se tiene una elevada densidad de fallas geológicas activas (Brown, Jr Ward, & Plafker, 1973), según Segura, Bungum, Lindholm, and Hernández (1999) las fallas sísmicas locales en términos estadísticos, generan el 59% de la amenaza sísmica total en Managua. El 41% restante resulta de la zona de subducción, de otras zonas en la cadena volcánica y de la zona montañosa de Nicaragua. 9

Según Schmoll, Knishinsky, and Dobrovolny (1975) la importancia de consideraciones geológicas para la reconstrucción de Managua fue absolutamente necesaria luego del terremoto de 1972. Como acción inmediata, las autoridades competentes de ese entonces encargaron un mapa de fallas y de la amenaza sísmica, que fue presentado junto con la matriz de planeación, por Consultants (1975) la Vice Ministro de Planificación Urbana (Figura 3.4).

Figura 3.4. Mapa de zonas de fallas geológicas principales. (Consultants, 1975)

Es de considerable importancia que el Mapa de Fallas Geológicas de Managua se encuentre en constante actualización, una de estas se realizó en Abril del 2002, sin embargo en este estudio no se logró realizar un mapa geologico-estructural el cual exige parámetros geométricos y estructurales (desplazamientos, grado de actividad, buzamiento, etc.) de las fallas o su edad relativa, se limitaron a presentar la ubicación y trayectoria de las fallas. (figura 3.5) (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002)

Figura 3.5. Zonas de fallas y áreas blancas que requieren futura investigación. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002)

10

3.7.

FALLAS SÍSMICAS Y AMENAZA VOLCÁNICA

W. Strauch (2000) detalla que es requerido pensar en la posibilidad de una futura actividad volcánica en las fallas sísmica principal, ubicadas en el mismo centro de Managua. El cráter Tiscapa es un ejemplo de la ocurrencia de un centro volcánico en una falla sísmica activa (Figura 3.6). Se afirma que la interrelación entre actividad volcánicas y tectónica-sísmica en la cadena volcánica de Nicaragua, fue demostrada claramente durante la erupción del volcán Cerro negro, en Agosto de 1999, cuando se produjeron sismos destructivos de magnitudes hasta 5 Richter (I. Strauch, Tenorio, & Bodán, 1999). Procesos similares podrían ocurrir en Managua.

Figura 3.6. Fotografía histórica de la Laguna de Tiscapa. (W. Strauch, 2000)

3.8.

HISTORIA SÍSMICA Y SISMICIDAD ACTUAL

3.8.1. Los Terremotos de 1931 y de 1972 El terremoto de 1931 revelo la existencia de una falla geológica que pasa por el estadio Nacional, sin embargo esto no estimuló la realización de muchos estudios técnicos. La fuente científica más importante sobre este evento corresponde a Sultán (1931). El terremoto de 1972, uno de los desastres más importantes acaecidos después de la adopción de la nueva teoría general de las placas tectónicas, alentó una gran cantidad de estudios a nivel nacional e internacional. El centro de Managua fue completamente destruido por un sismo de magnitud 5.6 Richter (6.2 MS), en la figura 3.7 se muestran algunas fotos representativas.

11

Figura 3.7. Fotos del terremoto de Managua, 1972 colección Steinbrugge. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002)

El acelerógrafo instalado en la refinería ESSO midió un nivel de aceleración horizontal de 0.39 g (1g corresponde a la gravedad de la tierra) a una distancia de 4 km de la falla Tiscapa y, aproximadamente, 7 km del hipocentro (Figura 3.8)

Figura 3.8. Acelerograma del terremoto de 1972. (W. Strauch, 2000)

12

Debe asumirse que la aceleración en la cercanía inmediata de las fallas activadas fue mucho más alta y podría haber alcanzado o sobrepasado 1g, en el centro de Managua. Los registros análogos fueron digitalizados, y se conservan en la base de datos de INETER. En el acelerógrafo de la UNAN, la aceleración máxima fue de 0.6, registrada durante una réplica fuerte que ocurrió en marzo de 1973. Esta amplitud tan alta se explica por la cercanía del epicentro del sismo al sitio de la falla activada. (Figura 3.9)

Figura 3.9. Epicentros de las réplicas del terremoto de Managua de 1972. (W. Strauch, 2000)

3.9.

SISMICIDAD ACTUAL

La Red Sísmica Nacional fue establecida en el año 1975, y ha localizado hasta ahora (con interrupción de 1985 a 1991), cerca de 25 mil sismos en Nicaragua, o sea, más de 1500 sismos por año (Figura 3.10). INETER desarrolla y mantiene esta red sísmica de 27 estaciones telemétricas. El registro digital funciona en la Central Sísmica, en Managua. Además existe una red de 20 estaciones sísmicas acelerográficas digitales ubicadas en Managua y en las cabeceras departamentales.

13

Figura 3.10. Red Sísmica Nacional de Nicaragua. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002)

La mayoría de los epicentros están localizados en el Océano Pacifico, donde chocan las placas tectónicas Coco y Caribe, En la figura 3.11 se aprecian los sismos superficiales (color rojo-hasta 40 km) en esta zona ocurren mar adentro. Los sismos de profundidad intermedia (amarillo) y muy profundos (azul-hasta 250 km) se dan más cerca de la costa del océano o directamente debajo de ella.

Figura 3.11. Vista Tridimensional de la sismicidad de Nicaragua. (Instituto Nicaragüense de Estudios Territoriales (INETER), 2002)

14

4. REVISIÓN DE LITERATURA

4.1. GENERALIDADES En este capítulo se abordaran algunos conceptos fundamentales relacionados al origen de los sismos. 4.1.1. Causas de los sismos Los fenómenos causantes de que la tierra tiemble son variados, en dependencia de esto se identifican tres clases: los sismos de origen volcánico, los artificialmente producidos por el hombre y los de origen tectónico. Este último es el más devastador y es por ello su mayor interés en la ingeniería (Goytia & Villanueva, 2001). 4.1.2. Tectónica de Placas La explicación del origen de los sismos puede brindarse recurriendo a la teoría de la tectónica de placas; esta sugiere que la corteza terrestre, la litosfera, se encuentra fraccionada en seis placas continentales (África, América, Antártida, Australia, Europa y la placa del Pacifico) y alrededor de 14 placas subcontinentales (Caribe, Coco, Nazca, Filipinas, etc.) (Figura 4.1). La teoría de la tectónica de placas fue desarrollada basándose en datos instrumentales hacia finales de la década de 1950. Estos datos mostraron que la mayor ocurrencia de terremotos se da en zonas estrechas y definidas, sugiriendo que la mayoría de sismos registrados resultan del choque entre estas placas. La explicación de la causa de este choque es el equilibrio térmico de los materiales que componen la Tierra, la explicación básica es la siguiente: El material caliente (en las profundidades) posee una densidad menor al material relativamente frio (cerca de la corteza terrestre), lo que hace que tienda a subir, mientras que el material de la superficie una vez frio tiende a bajar por acción de la gravedad. Este proceso cíclico es denominado convección. Las corrientes convectivas generan esfuerzos de corte en la base de las placas, provocando su movimiento en todas direcciones.(ROSENBLUETH, 1992)

15

Figura 4.1. Principales zonas tectónicas, lomos oceánicos y zonas de subducción. (Kramer, 1996)

4.1.3. Fallas Geológicas Las fuentes generadoras de sismos generalmente esta asociadas a anomalías geológicas capaces de acumular gran cantidad de energía producto de la naturaleza dinámica de la corteza terrestre. Estas estructuras son las fallas geológicas las cuales se manifiestan como fracturas ocasionadas por el desplazamiento relativo de los dos lados de la ruptura. La longitud de estas fallas puede ser de varios metros hasta cientos de kilómetros y ser extendidos desde la superficie a varias decenas de kilómetros de profundidad.

Figura 4.2. Falla de San Andrés (falla por desgarramiento). (ROSENBLUETH, 1992)

La presencia de estas fallas en la superficie no implica necesariamente que el área posee actividad sísmica, así como la inexistencia de las mismas no indica que el área sea asísmica, por cuanto muchas veces estas fracturas no llegan a aflorar en la superficie. 16

4.1.4. Tipos de fallas En dependencia de su movimiento se dan tres casos: normal, inversa y de desgarradura (Figura 4.3), las fallas normales son propias de la zona de tracción las fallas inversas corresponden a zonas de compresión y las fallas por desgarramiento implican grandes desplazamientos laterales entre d os placas en contacto, la falla de San Andrés es un ejemplo ilustrativo de este tipo (Figura 4.2) (a)

(b)

(c)

Figura 4.3. Tipos de fallas según su patrón de ruptura: (a) Fallas normales, (b) Fallas inversas, (c) Fallas transformantes (Kramer, 1996).

4.2. ONDAS SÍSMICAS La liberación de energía en el foco o hipocentro del sismo, se propaga en forma de ondas elásticas de deformación. Se tiene la primicia que las deformaciones generadas por el paso de una onda son elásticas, en base a esto las velocidades de propagación son determinadas mediante los parámetros de modulo elástico y densidad de los materiales a través de los cuales viaja la onda. Las ondas sísmicas son clasificadas según su naturaleza en ondas de cuerpo y ondas de superficie. (Bolt, 1993) 4.2.1. Ondas de cuerpo Las ondas de cuerpo han sido de gran utilidad para propósitos de estudios de la estructura interna de la tierra y en la actualidad tienen un gran número de aplicaciones en ingeniería las cuales se dividen en ondas P y Ondas S. Las ondas P, denominadas también primarias, longitudinales, de compresión o dilatación; producen movimientos de partículas en la misma dirección de la propagación, alternando compresión y dilatación del medio (Figura 4.4). (Bolt, 1993) 17

Figura 4.4. Deformaciones producidas por las ondas P (Bolt, 1993)

Las ondas S, también denominadas ondas secundarias, transversales y de corte; producen un movimiento de partículas en sentido perpendicular a la dirección de la propagación (Figura 4.5).

Figura 4.5. Deformaciones producidas por las ondas S (Bolt, 1993)

Generalmente cuando se da un sismo, las ondas P se registran en primera instancia, seguido llegan las ondas S. Las ondas P se propagan a través de medios sólidos y líquidos, mientras que las ondas S, se propagan exclusivamente en medios solidos por cuanto los líquidos no presentan resistencia al corte. 4.2.2. Ondas Superficiales Estas son denominadas de esta manera por cuanto su movimiento es restringido a las cercanías de la superficie terrestre. Se clasifican en dos tipos: ondas Love y ondas Rayleigh. El movimiento de las ondas L (Figura 4.6), es similar al de las ondas S que no poseen componente vertical esta mueve la superficie del suelo de lado a lado sobre un plano horizontal y en sentido perpendicular a la dirección de propagación.

18

Figura 4.6. Deformaciones producidas por las ondas L (Bolt, 1993)

El movimiento de las partículas en las ondas R es elíptico y tiene lugar en planos perpendiculares a la superficie libre. (Figura 4.7)

Figura 4.7. Deformaciones producidas por las ondas R (Bolt, 1993)

En general, las ondas Love son más veloces que las ondas Rayleigh, pero ambas se propagan a menor velocidad que las ondas de cuerpo. El intervalo de llegada entre las ondas puede observarse en algunos acelerogramas tal es el caso del acelerograma del terremoto en Kermadec (Figura 4.8).

Figura 4.8. Terremoto de Kermadec de 11 de Junio de 1957 (Kramer, 1996)

19

4.3. RESPUESTA SÍSMICA EN 1D Las simulaciones de respuesta sísmica del suelo frente a un movimiento sísmico han sido realizadas mediante programas informáticos que utilizan varias hipótesis simplificadoras. Uno de estos primeros programas fue SHAKE (Schnabel, 1972) este programa analiza la respuesta de un sistema roca-suelo teniendo como premisa que los estratos son horizontalmente infinitos por tanto no se considera la variación lateral del suelo (heterogeneidad) y estos sometidos a una incidencia vertical de ondas sísmicas S. También tiene como supuesto que el comportamiento cíclico del suelo puede aproximarse mediante un modelo lineal equivalente. Está basado en las soluciones de propagación de onda de Kanai (1951), Roesset and Whitman (1969) y Tsai and Housner (1970). Otros programas desarrollados que aplican el método lineal equivalente son EERA (Bardet, Ichii, & Lin, 2000) y DEEPSOIL (Hashash, Groholski, Phillips, & Park, 2009). Sin embargo si se desea considerar el efecto de heterogeneidad de suelos serán requeridos algoritmos más avanzados tales como elemento finito. El modelo 1D equivalente representa la respuesta tensión deformación del suelo basándose en el modelo Kelvin-Voigt (figura 4.9).

Figura 4.9. Representación esquemática del modelo tensión-deformación utilizada en el modelo lineal equivalente. (J.F. Semblat & Pecker, 2009).

La tensión de corte 𝜏 está en dependencia de la deformación de corte 𝛾 y de su primera derivada mediante la ecuación: 𝜏 = 𝐺𝛾 + 𝜂𝛾̇

(4.1)

El análisis unidimensional lineal equivalente de la respuesta de suelo es representado en la figura 4.10. El depósito de suelo está formado por N estratos de suelo, horizontalmente infinitos, de espesor h, y con unas propiedades características de cada estrato: densidad, módulo de corte y razón de amortiguamiento.

20

Figura 4.10. Depósito de suelo estratificado para su estudio mediante el método unidimensional. (Schnabel, 1972)

La ecuación unidimensional del movimiento provocado por la propagación de ondas de corte es: 𝜕2 𝑢

𝜕𝜏

𝜌 𝜕𝑡 2 = 𝜕𝑧

(4.2)

Si se supone que cada estrato se comporta como un sólido Kelvin-Voigt, la ecuación anterior usando las ecuaciones, se transforma en: 𝜕2 𝑢

𝜕2 𝑢

𝜕3𝑢

𝜌 𝜕𝑡 2 = 𝐺 𝜕𝑧 2 + 𝜂 𝜕𝑧 2𝜕𝑡

(4.3)

Expresando la ecuación en función de la expresión del desplazamiento de un movimiento armónico se obtiene: (𝐺 + 𝑖𝜔𝜂)

𝑑2 𝑢

= 𝜌𝜔2 𝑈

𝑑𝑧 2

(4.4)

Esta ecuación admite la siguiente ecuación general: 𝑈(𝑥) = 𝐸𝑒 𝑖𝑘

∗𝑧

+ 𝐹𝑒 −𝑖𝑘

∗𝑧

(4.5)

Donde E, F son las amplitudes del desplazamiento y 𝑘 ∗ es el número de onda complejo cuya expresión es: 𝜌𝜔 2

𝑘 ∗ = 𝐺+𝑖𝜔𝜂 =

𝜌𝜔 2

(4.6)

𝐺∗

La solución de la ecuación es: 𝑢(𝑧, 𝑡) = (𝐸𝑒 𝑖𝑘

∗𝑧

∗

+ 𝐹𝑒 −𝑖𝑘 𝑧 )𝑒 𝑖𝜔𝑡

(4.7) 21

Y la tensión de corte correspondiente es: ∗

∗

𝜏(𝑧, 𝑡) = 𝑖𝑘 ∗ 𝐺 ∗ (𝐸𝑒 𝑖𝑘 𝑧 − 𝐹𝑒 −𝑖𝑘 𝑧 )𝑒 𝑖𝜔𝑡

(4.8)

Los desplazamientos en la base (𝑧 = 0) y en el techo (𝑧 = ℎ𝑚 ) del estrato m son: 𝑢𝑚 (0, 𝑡) = (𝐸𝑚 + 𝐹𝑚 )𝑒 𝑖𝜔𝑡

𝑢(𝑧, 𝑡) = (𝐸𝑚 𝑒 𝑖𝑘

∗ℎ 𝑚

+ 𝐹𝑚 𝑒 −𝑖𝑘

∗ℎ 𝑚

)𝑒 𝑖𝜔𝑡

(4.9)

Donde 𝐸𝑚 y 𝐹𝑚 son las amplitudes del desplazamiento para un determinado estrato m. Las tensiones de corte en base y techo del estrato m son respectivamente: ∗ ∗ (𝐸 𝜏𝑚 (0, 𝑡) = 𝑖𝑘𝑚 𝐺𝑚 𝑚 + 𝐹𝑚 )𝑒 𝑖𝜔𝑡 ∗ 𝐹𝑚 𝑒 −𝑖𝑘𝑚 ℎ𝑚 )𝑒 𝑖𝜔𝑡

∗ ∗ 𝜏𝑚 (ℎ𝑚 , 𝑡) = 𝑖𝑘𝑚 𝐺𝑚 (𝐸𝑚 𝑒 𝑖𝑘

∗ℎ 𝑚

− (4.10)

En el contacto entre dos estratos sucesivos 𝑚 y 𝑚 + 1 los desplazamientos y las tensiones de corte deben ser también continuos. Esto implica que: 𝑢𝑚 (ℎ𝑚 , 𝑡) = 𝑢𝑚+1 (0, 𝑡)

𝜏𝑚 (ℎ𝑚 , 𝑡) = 𝜏𝑚+1 (0, 𝑡)

(4.11)

Usando las ecuaciones se pueden relacionar los coeficientes 𝐸𝑚 y 𝐹𝑚 : 𝐸𝑚+1 + 𝐹𝑚+1 = 𝐸𝑚 𝑒 𝑖𝑘 𝐸𝑚+1 − 𝐹𝑚+1 = 𝑘 ∗

∗ℎ 𝑚

∗ 𝐺∗ 𝑘𝑚 𝑚

∗ 𝑚+1 𝐺𝑚+1

+ 𝐹𝑚 𝑒 −𝑖𝑘𝑚

(𝐸𝑚 𝑒 𝑖𝑘

∗ℎ 𝑚

∗

ℎ𝑚

(4.12)

− 𝐹𝑚 𝑒 −𝑖𝑘𝑚

∗

ℎ𝑚

)

(4.13)

De las ecuaciones se obtienen los siguientes algoritmos que permiten obtener las amplitudes del desplazamiento de los estratos superiores, 𝐸𝑚+1 y 𝐹𝑚+1 en función de las amplitudes de los estratos inferiores, 𝐸𝑚 y 𝐹𝑚 : 1

∗

1

∗

∗ )𝑒 𝑖𝑘𝑚 𝐸𝑚+1 = 2 𝐸𝑚 (1 + 𝛼𝑚 ∗ )𝑒 𝑖𝑘𝑚 𝐹𝑚+1 = 2 𝐸𝑚 (1 − 𝛼𝑚

1

∗

1

∗

ℎ𝑚

∗ )𝑒 −𝑖𝑘𝑚 + 2 𝐸𝑚 (1 − 𝛼𝑚

ℎ𝑚

∗ )𝑒 −𝑖𝑘𝑚 + 2 𝐸𝑚 (1 − 𝛼𝑚

ℎ𝑚

(4.14)

ℎ𝑚

(4.15)

∗ Donde 𝛼𝑚 es la razón de impedancia en forma compleja en la interfaz de dos estratos 𝑚 y 𝑚 + 1 ∶

∗ 𝛼𝑚 = 𝑘∗

∗ 𝐺∗ 𝑘𝑚 𝑚

∗ 𝑚+1 𝐺𝑚+1

= √𝜌

∗ 𝜌𝑚 𝐺𝑚

(4.16)

∗ 𝑚+1 𝐺𝑚+1

22

El algoritmo anterior se inicia en la superficie libre donde no se presenta tensión de corte, es decir: 𝜏1 (0, 𝑡) = 𝑖𝑘1∗ 𝐺1∗ (𝐸1 − 𝐹1 )𝑒 𝑖𝜔𝑡 = 0

(4.17)

Lo que implica que: 𝐸1 = 𝐹1

(4.18)

El algoritmo obtenido se aplica entonces sucesivamente de los estratos 2 a 𝑚. La función de transferencia que relaciona los desplazamientos de los estratos 𝑚 y 𝑛 se define como: 𝐴𝑚𝑛 (𝜔) =

𝑢𝑚 𝑢𝑛

=

𝐸𝑚 +𝐹𝑚

(4.19)

𝐸𝑛 +𝐹𝑛

La velocidad y la aceleración se relacionan con el desplazamiento a través de: 𝑢̇ (𝑧, 𝑡) =

∂u 𝜕𝑡

= 𝑖𝜔𝑢(𝑧, 𝑡) 𝑦 𝑢̈ (𝑧, 𝑡) =

𝜕2 𝑢 𝜕𝑡 2

= −𝜔2 𝑢(𝑧, 𝑡)

(4.20)

Por tanto la función de transferencia también puede ser expresada en función de las velocidades y las aceleraciones en los estratos 𝑚 y 𝑛: 𝐴𝑚𝑛 (𝜔) =

𝑢𝑚 𝑢𝑛

=

𝑢̇ 𝑚 𝑢̇ 𝑛

=

𝑢̈ 𝑚 𝑢̈ 𝑛

=

𝐸𝑚 +𝐹𝑚

(4.21)

𝐸𝑛 +𝐹𝑛

Por lo tanto la función de transferencia compara el desplazamiento, velocidad o aceleración entre base y techo de un estrato. La deformación de corte a una profundidad 𝑧 y tiempo t puede ser obtenida mediante la ecuación: 𝛾(𝑧, 𝑡) =

𝜕𝑢 𝜕𝑧

= 𝑖𝑘 ∗ (𝐸𝑒 𝑖𝑘

∗𝑧

∗

− 𝐹𝑒 −𝑖𝑘 𝑧 )𝑒 𝑖𝜔𝑡

(4.22)

Y la tensión de corte correspondiente es: 𝜏(𝑧, 𝑡) = 𝐺 ∗ 𝛾(𝑧, 𝑡)

(4.23)

En los análisis de respuesta de sitio es requerido definir la localización y tipo de movimiento sísmico, distinguiéndose los siguientes términos: movimiento en la superficie libre de un depósito de suelo (free surface motion), movimiento en el sustrato rocoso (bedrock motion) y movimiento en afloramiento de roca (rock outcropping motion) (Figura 4.11).

23

Figura 4.11. Terminología utilizada en análisis de respuesta de sitio. (J.F. Semblat & Pecker, 2009)

Como se muestra en la figura, la onda incidente S (que se propaga verticalmente hacia arriba a través del sustrato rocoso) tiene una amplitud 𝐸𝑁 . El movimiento, ya en el sustrato rocoso tiene una amplitud 𝐸𝑁 + 𝐹𝑁 , en el techo del sustrato rocoso y justo debajo de los estratos del suelo. Debido a que en la superficie libre no se presentan tensiones de corte, sustituyendo la ecuación en la ecuación se obtiene que la amplitud del movimiento en un afloramiento es: 2𝐸𝑁 . Por lo tanto la función de transferencia relativa a los movimientos en un afloramiento respecto al sustrato rocoso es: 𝐴𝑀𝑁 (𝜔) = 𝐸

2𝐸𝑁

(4.24)

𝑁 +𝐹𝑁

4.3.1. La aproximación lineal equivalente La aproximación lineal equivalente consiste en modificar el modelo Kelvin-Voigt para tener en cuenta algunos tipos de comportamiento no lineales del suelo. El módulo de corte lineal equivalente es tomado igual al módulo de corte secante el cual es dependiente de la amplitud de la deformación de corte 𝛾𝑐 según: 𝜏

𝐺𝑠𝑒𝑐 = 𝛾𝑐

(4.25)

𝑐

En los análisis de respuesta de sitio el comportamiento del material se especifica al estudiar las variaciones del módulo de corte secante y de la razón de amortiguamiento respecto a la deformación de corte: Las curvas 𝐺𝑠𝑒𝑐 − 𝛾 no pueden tener formas arbitrarias ya que derivan de curvas 𝜏 − 𝛾. Esto implica que cualquier restricción en las relaciones 𝜏 − 𝛾 debe ser reflejada también en restricciones en las relaciones 𝐺𝑠𝑒𝑐 − 𝛾. Por ejemplo, existe el fenómeno del reblandecimiento en los suelos, es decir, la disminución de corte con la deformación. El estudio de este tipo de fenómeno requiere de técnicas numéricas especiales y complicadas por cuanto su estudio con las técnicas de trabajo habituales en la respuesta de sitio da lugar a

24

efectos numéricos tales como que la solución depende fuertemente de la discretización espacial. La exclusión del reblandecimiento implica que: 𝑑𝜏

= 𝐺𝑠𝑒𝑐 (𝛾) + 𝑑𝛾

𝑑𝐺𝑠𝑒𝑐 𝑑𝛾

𝛾 ≥ 0 lo que se traduce en

∆𝐺𝑠𝑒𝑐 𝐺𝑚𝑎𝑥

≥−

𝐺𝑠𝑒𝑐 (𝛾) ∆𝛾 𝐺𝑚𝑎𝑥

𝛾

(4.26)

Donde ∆𝐺𝑠 es la disminución en el valor de 𝐺𝑠 correspondiente a un aumento ∆𝛾 en la deformación y 𝐺𝑚𝑎𝑥 es el valor máximo del módulo de corte máximo. En el modelo lineal equivalente la hipótesis inicial es la dependencia del módulo de corte y la razón de amortiguamiento con la deformación de corte. Para obtener las dos primeras se realizan iteraciones hasta que el valor obtenido sea consistente con el nivel de deformación inducido en cada estrato. Una breve descripción es la siguiente: Se parte de los valores iniciales de módulo de corte 𝐺0 y razón de amortiguamiento 𝜉0 correspondientes a niveles de deformación pequeña, luego es calculada la deformación de corte efectiva para cada estrato, según: 𝛾𝑒𝑓𝑓 = 𝑅𝛾 𝛾𝑚𝑎𝑥

(4.27)

Donde 𝑅𝛾 es una razón que relaciona la deformación de corte efectiva con la deformación de corte máxima y dependiente de la magnitud del terremoto según la ecuación. A partir del valor 𝛾𝑒𝑓𝑓 se calculan los nuevos valores 𝐺 𝑖+1 y 𝜉 𝑖+1 con una nueva iteración. Finalmente se habrán de repetir los pasos anteriores hasta alcanzar la convergencia (Figura 4.12).

Figura 4.12. Proceso de iteración para el valor del módulo de corte y la razón de amortiguamiento en el análisis equivalente lineal. (J.F. Semblat & Pecker, 2009)

25

4.4. RESPUESTA SÍSMICA EN 2D 4.4.1. Ecuación de Equilibrio Salençon (2001), establece que la ecuación de equilibrio es obtenida por la expresión que la suma de la razón de trabajo por fuerzas externas y la razón de trabajo por fuerzas internas es igual a la razón de trabajo por cantidades de aceleración, es decir: ̂ )𝑑Ω + ∫ 𝑓. 𝑈 ̂ 𝑑Ω + ∫ 𝑇 𝑑 . 𝑈 ̂ 𝑑𝒶 = ∫ 𝜌𝒶 . 𝑈 ̂ 𝑑Ω − ∫Ω 𝜎: 𝑑̂(𝑈 Ω δΩ Ω

(4.28)

Tomando en cuenta la simetría de la fuerza de tensión, el primer término se convierte en: ̂ )𝑑Ω = ∫ 𝜎: ∇𝑈 ̂ 𝑑Ω ∫Ω 𝜎: 𝑑̂ (𝑈 Ω

(4.29)

Y la integración por partes: ̂ )𝑑Ω = ∫ (∇. 𝜎) . 𝑈 ̂ 𝑑Ω − ∫ (𝜎. 𝑛) . 𝑈 ̂ 𝑑𝒶 − ∫Ω 𝜎: 𝑑̂(𝑈 Ω ∂Ω

(4.30)

Introduciendo la Ecuación 4.30 en la ecuación 4.28 y escribiendo la ecuación para ̂ , obtenemos la ecuación de equilibrio: cada campo de velocidad virtual 𝑈 ∇. 𝜎 + 𝑓 = 𝜌𝒶 𝜎. 𝑛 = 𝑇 𝑑

𝑒𝑛 Ω

(4.31)

𝑒𝑛 ∂Ω

(4.32)

Donde 𝜎 es la fuerza de tensión Cauchy, 𝑓 la fuerza del cuerpo y 𝑇 𝑑 la fuerza de superficie. 4.4.2. Ecuación Constitutiva La ley constitutiva elástico lineal isotrópico es: 𝜎 = 𝜆 (tr𝜀) Π + 2𝜇𝜀

(4.33)

Donde 𝜆 y 𝜇 representan la constante de Lamé, Π la fuerza unitaria de segundo orden (i.e son las componentes de 𝜀𝑖𝑗 ), y 𝜀 el tensor de deformación relacionado al campo de desplazamiento 𝑢 por: 1

𝜀 = 2 ( 𝑡∇𝑢 + ∇𝑢)

(4.34) 26

Es decir, en término de las componentes: 𝜕𝑢

1 𝜕𝑢

𝜀𝑖𝑗 = 2 (𝜕𝑥 𝑖 + 𝜕𝑥𝑗) 𝑗

(4.35)

𝑖

En la ecuación (4.33), los términos tr𝜀 corresponden a la deformación volumétrica: 𝜕𝑢

tr𝜀 = ∑3𝑖=1 𝜀𝑖𝑗 = ∑3𝑖=1 𝜕𝑥𝑖

(4.36)

𝑖

La constante de Lamé, 𝜆 y 𝜇, se definen completamente con el comportamiento elástico lineal isotrópico de los materiales. En la tabla 4.1, recordamos la relación entre 𝜆,𝜇, del módulo de Young´s 𝐸, la relación de Poisson 𝑣 y el módulo de compresibilidad 𝐾. Tabla 4.1. Relaciones entre los parámetros definidos por el comportamiento elástico lineal isotrópico de los materiales. (J.F. Semblat & Pecker, 2009)

𝑬 𝐸

𝝂 𝜈

𝜇(3𝜆 + 2𝜇) 𝜆+𝜇

𝜆 2(𝜆 + 𝜇)

𝑬, 𝝂 𝜆, 𝝁

𝝀 𝜈𝐸 (1 + 𝜈)(1 − 2𝜈) 𝜆

𝝁=𝑮 𝐸 2(1 + 𝜈) 𝜇

𝑲 𝐸 3(1 − 2𝜈) 3𝜆 + 2𝜇 3

4.4.3. Ecuación de Equilibrio en términos de desplazamiento Incluyendo la ley constitutiva (Ecuación 4.33) en la ecuación de equilibrio (4.31) y tomando en cuenta la definición de tensor de deformaciones (Ecuación. 4.34), obtenemos la ecuación de equilibrio expresada en función al campo de desplazamiento 𝑢: (𝜆 + 𝜇)[∇(∇. 𝑢)] + 𝜇∆𝑢 = 𝜌

𝜕2 𝑢 𝜕𝑡 2

−𝑓

(4.37)

Donde los operadores ∇ y ∆ representan el gradiente y Laplaciana respectivamente. En términos de los componentes, el rendimiento de la ecuación (4.37): (𝜆 + 𝜇)

𝜕 𝜕𝑥𝑖

𝜕𝑢

𝜕2 𝑢

(∑3𝑗=1 𝜕𝑥𝑗 ) + 𝜇 ∑3𝑗=1 𝜕𝑥 2 𝑗 = 𝜌 𝑗

𝑗

𝜕2 𝑢𝑖 𝜕𝑡 2

− 𝑓 𝑖 ; 𝑖 = 1,3

(4.38)

4.4.4. Descomposición del campo de desplazamiento Considerando el teorema de Poisson, siempre es posible descomponer un campo vectorial 𝑢 como la suma de un gradiente de un potencial escalar ø y el rotacional de un potencial vectorial 𝜓: 27

𝑢 = ∇(𝜙) + ∇ Λ 𝜓

(4.39)

Es decir, en términos de componentes: 𝜕∅

𝜕𝜓𝑘

𝑖

𝜕𝑥𝑗

𝑢𝑖 = 𝜕𝑥 +

𝜕𝜓

− 𝜕𝑥 𝑗 , para la permutación circular 𝑖, 𝑗, 𝑘 𝑘

(4.40)

Examinando la ecuación (4.39) y a partir de los tres componentes del vector 𝑢, introducimos cuatro componentes por la descomposición en potenciales: el potencial escalar y los tres componentes del potencial vectorial. Es entonces necesario especificar una condición adicional. Tal condición se obtiene generalmente por escrito, que la divergencia del potencial vectorial 𝜓 es cero: ∇. 𝜓 = 0

(4.41)

Esto significa que el potencial vectorial 𝜓 no llevó a ningún cambio de volumen. 4.4.5. Ecuaciones de ondas desacopladas Incluyendo la ecuación 4.39 en la ecuación de equilibrio (Ecuación 4.37) y tomando en cuenta la ecuación 4.41 obtenemos las siguientes ecuaciones para los potenciales: 1 𝜕2 ∅

∆(𝜙) = 𝑉 2 𝑃

𝜕𝑡 2

−𝐹

2 1 𝜕 𝜓

∆ (𝜓) = 𝑉 2 𝑆

𝜕𝑡 2

(4.42)

−𝐺

(4.43)

En el que se consideran las notaciones: 𝜆+2𝜇

𝑉𝑃 = √

𝜇

Y 𝑉𝑆 = √𝜌

𝜌

(4.44)

Y aplicando el teorema de Poisson de 𝑓: 𝑓 = 𝑉𝑃2 ∇(𝐹) + 𝑉𝑆2 ∇ Λ (𝐹)

(4.45)

En términos de componentes, las ecuaciones (4.42) y (4.43) puede escribirse: ∑3𝑗=1 ∑3𝑗=1

𝜕2 𝜙 𝜕𝑥 2

𝑗

𝜕2 𝜓𝑖 𝜕𝑥 2

𝑗

1 𝜕2 𝜙

= 𝑉2 𝑃

𝜕𝑡 2

−𝐹

1 𝜕2 𝜓𝑖

= 𝑉2 𝑆

𝜕𝑡 2

(4.46)

− 𝐺𝑖 , 𝑖 = 1,3

(4.47)

28

4.4.6. Ondas de cuerpo Omitiendo la condición límite, la ecuación (4.42) y (4.43) rigen la propagación de dos diferentes velocidades de propagación de ondas 𝑉𝑃 y 𝑉𝑆 . Estas ondas son llamadas ondas de cuerpo y son de dos tipos: Ecuación (4.42) rige la propagación de la onda de presión que también se llama onda-P, Ecuación (4.43) rige la propagación de la onda de corte que también se llama ondaS. La polarización de estas ondas es ilustrada en la siguiente sección para ondas planas. La comparación de las ecuaciones (4.44) muestra que la velocidad 𝑉𝑃 de la onda P es mayor que la velocidad 𝑉𝑆 de la onda S, además J.F. Semblat and Pecker (2009) definen valores típicos para las velocidades 𝑉𝑃 y 𝑉𝑠 según el tipo de material (Tabla 4.2). La relación 𝑉𝑃 ⁄𝑉𝑆 puede ser fácilmente expresada como una función de la relación de Poisson: 𝑉𝑃 𝑉𝑆

2−2𝜈

= √1−2𝜈

(4.48)

Figura 4.13. Relación entre la velocidad de la onda P y la velocidad de la onda S.

La figura 4.13 muestra las variaciones de la relación entre la velocidad de la onda P y la velocidad de la onda S. La velocidad de la onda P es siempre mayor que la velocidad de la onda S y la proporción aumenta rápidamente para relaciones de Poisson por encima de 0.4.

29

Tabla 4.2. Valores típicos de velocidades para las ondas P y S para diferentes tipos de suelo. (J.F. Semblat & Pecker, 2009)

Material Arcilla Arena arriba del nivel freático Arena debajo del nivel freático Marga Roca

𝑽𝒔 (𝒎⁄𝒔) 100-200 200-400 200-400 400-600 >800

𝑽𝒑 (𝒎⁄𝒔) 1500 400-800 1500-1700 1500-2000 >2000

4.4.7. Propagación de Ondas en medio anisotrópico En el caso de la elasticidad lineal en pequeñas deformaciones, la ley constitutiva tiene la siguiente forma: 𝜎 = 𝐶: 𝜀

(4.49)

Donde 𝜎 es la fuerza de tensión, 𝐶 es el tensor de cuarto orden y 𝜀 es el tensor de deformaciones. En el caso isotrópico, el tensor de elasticidad conduce a la misma expresión como en la ecuación (4.33), ya que se puede expresar como: 𝐶𝑖𝑗𝑘𝑙 = 𝜆𝛿𝑖𝑗 𝛿𝑘𝑙 + 𝜇(𝛿𝑖𝑘 𝛿𝑗𝑙 + 𝛿𝑖𝑙 𝛿𝑗𝑘 )

(4.50)

Para un medio anisotrópico, 𝐶 puede tener hasta 21 componentes independientes. Medios geológicos en ocasiones son considerados como medios ortotrópicos (también llamados medios transversalmente anisotrópico). Suponiendo el eje ortotrópico como 𝑒𝑧 , todas las direcciones de los planos perpendiculares a 𝑒𝑧 son entonces equivalente. Entre los 21 componentes del tensor de elasticidad, solo cinco son independientes: 𝐶𝑥𝑥𝑥𝑥 = 𝐶𝑦𝑦𝑦𝑦 , 𝐶𝑦𝑦𝑦𝑦 , 𝐶𝑥𝑥𝑦𝑦 , 𝐶𝑦𝑦𝑧𝑧 = 𝐶𝑥𝑥𝑧𝑧 𝑦 𝐶𝑥𝑧𝑥𝑧 = 𝐶𝑦𝑧𝑦𝑧 (Bourbié, Coussy, Zinszner, Junger, & Miguel, 1987) todos los demás componentes son cero excepto 𝐶𝑥𝑦𝑥𝑦 que está escrito: 1

𝐶𝑥𝑦𝑥𝑦 = 2 (𝐶𝑥𝑥𝑥𝑥 − 𝐶𝑥𝑥𝑦𝑦 )

(4.51)

Finalmente obtenemos cuatro valores de velocidad independientes (Bourbié, et al., 1987)

30

Para ondas P: 𝐶𝑧𝑧𝑧𝑧

𝑉𝑃 = √ 𝐶𝑥𝑥𝑥𝑥

𝑉 =√ { 𝑃

𝜌

𝜌

(4.52)

𝐶𝑦𝑦𝑦𝑦

=√

𝜌

Para ondas S: 𝐶𝑦𝑧𝑦𝑧

𝑉𝑆 = √

𝐶𝑥𝑧𝑥𝑧

=√

𝜌

𝜌

(4.53)

𝐶𝑥𝑦𝑥𝑦

{

𝑉𝑆 = √

𝜌

En el caso de los medios anisotrópicos, la velocidad de la onda depende por lo tanto de la dirección de propagación. Su caracterización se puede realizar mediante el uso de mediciones de campo (Aki & Richards, 2002) o las pruebas de laboratorio (Nguyen, Semblat, Reiffsteck, & Lenti, 2008) 4.4.8. Comparación para casos de una y dos capas Con el fin de estimar la influencia de la estratificación del suelo en la amplificación de las ondas sísmicas, ahora se propone una comparación entre el caso de una sola capa y varios casos de dos capas. El caso de una sola capa (caso 1) corresponde a una sola capa de espesor h1=20 𝑚 y velocidad Vs1 = 200 𝑚/𝑠, que cubre un semi-espacio.

Figura 4.14. Función de transferencia a través de una capa elástica del suelo (caso 1) o a través de dos capas elásticas con diferentes velocidades (casos 2, 3 y 4). (J.F. Semblat & Pecker, 2009)

31

En los casos de dos capas (casos 2, 3 y 4), elegimos dos capas gruesas de 10 m correspondientes a la misma capa de espesor total que en el caso anterior. Las velocidades en cada capa se determinan a fin de aumentar progresivamente la relación de velocidad entre las dos capas y al mismo tiempo mantener la misma 𝑉 +𝑉 velocidad equivalente ( 1 2 2 = 200 𝑚⁄𝑠 como en el caso 1) Caso 2: 𝑉𝑆1 = 150 𝑚⁄𝑠 ; 𝑉𝑆2 = 250 𝑚⁄𝑠, Caso 3: 𝑉𝑆1 = 100 𝑚⁄𝑠 ; 𝑉𝑆2 = 300 𝑚⁄𝑠, Caso 4: 𝑉𝑆1 = 75 𝑚⁄𝑠 ; 𝑉𝑆2 = 325 𝑚⁄𝑠, Los resultados mostrados en la Figura. 4.14 indican que, para un espesor total de capa idéntica, las dos capas de los casos conducen a una amplificación más fuerte que en el caso de una sola capa. Por otra parte, cuanto mayor sea la relación de velocidad entre las dos capas, mayor es la amplificación. Por último, la frecuencia de la amplificación máxima disminuye a medida que aumenta la relación de velocidad. 4.4.9. Modelamiento de propagación de ondas Los problemas de propagación de ondas están caracterizados por varios fenómenos (Aki & Richards, 2002; Bourbié, et al., 1987) dispersión, esparcimiento, atenuación, conversión de onda, etc. Los parámetros que gobiernan todos estos fenómenos no siempre pueden ser estimados experimentalmente. A menudo es necesario para llevar a cabo modelos experimentales (modelo de materiales (Bodet et al., 2005), modelos a escala reducida (Coe, Prevost, & Scanlan, 1985; Chazelas, Abraham, & Semblat, 2001), etc.). También se puede considerar modelos numéricos y/o métodos inversos con el fin de caracterizar el material y las ondas que se propagan a través de él (Foti, 2000, 2003; Lai, Rix, Foti, & Roma, 2002). Los diferentes métodos numéricos que permiten el modelado de la propagación de ondas (diferencias finitas, elementos finitos, elementos de frontera, etc.) tienen diferentes ventajas y desventajas. El método de elementos finitos, por ejemplo, es conveniente tener una descripción confiable del comportamiento del material, pero no permite un fácil modelado de la propagación de ondas dentro del diámetro que tiene una extensión infinita. Por otra parte, la diferencia finita o los métodos de elemento finito inducen una dependencia de las características de propagación en algunos parámetros algorítmicos (esquema de integración en el tiempo, discretización espacial, etc.) (Bamberger, Chavent, & Lailly, 1980; Ihlenburg & Babuška, 1995; J.F Semblat & Brioist, 2000). Otros métodos numéricos (elementos de frontera, elementos espectrales) eliminan o limitan tales fenómenos (Beskos, 1997; M. Bonnet, 1999; Faccioli, Maggio, Paolucci, & Quarteroni, 1997; Komatitsch & Vilotte, 1998). Por lo tanto es necesario discutir las ventajas y desventajas de cada método numérico.

32

4.4.9.1. Dominio de Tiempo contra Dominio de Frecuencia La modelización numérica de ondas y vibraciones en los suelos se puede realizar ya sea por la integración dominio de tiempo o por métodos de solución de dominio de frecuencia. En el campo de la dinámica estructural, métodos modales, incluyendo superposición, se utilizan a menudo, pero no están plenamente adaptados para analizar la propagación de ondas (Clough & Penzien, 1993) Métodos de integración de tiempo permiten una implementación más fácil de leyes constitutivas dependientes del tiempo, así como las leyes constitutivas no lineales. Consisten en una estimación numérica de la solución en varias ocasiones a partir de las condiciones iniciales y las variaciones de tiempo de la carga. Estos métodos se basan en supuestos que permiten la aproximación de las derivadas temporales de los desplazamientos. Para leyes constitutivas lineales, la solución se obtiene mediante la resolución de un conjunto lineal de ecuaciones. Para problemas no lineales, en general es necesario implementar los procedimientos numéricos iterativos en cada intervalo de tiempo. Muchos algoritmos de integración de tiempo pueden ser considerados (Hughes, 1987), algunos de ellos se discutirán a continuación. Métodos de solución de dominio de frecuencia se utilizan ampliamente en los enfoques de los elementos de contorno (M. Bonnet, 1999). Son muy potentes en el campo de la dinámica / propagación de ondas desde la ecuación de movimiento se puede derivar fácilmente en el dominio de frecuencia (derivadas temporales expresan como productos simples por un factor de + -). 4.4.9.2. Señales reales o sintéticas Modelamiento numérico de los fenómenos de propagación de onda plantea la necesidad de la estimación de las propiedades dinámicas del medio, sino también las características de la carga. Diversos tipos de cargas pueden considerarse en función del objetivo del análisis (diseño, estudio paramétrico, etc):   

Señales sísmicas reales o vibraciones registrados en un sitio que se puede usar (directamente o después de deconvolución) como una carga / de entrada en un modelo numérico, Señales sintéticas calculadas a partir de los espectros estándar del sitio, Señales exclusivamente sintéticas determinadas a partir de expresiones analíticas simples.

4.4.9.3. Sismogramas Sintéticos En el caso de excitaciones sísmicas, el cálculo de acelerogramas reales por lo general requiere un conocimiento detallado de la fuente, la ruta y efectos de sitio. Algunos modelos empíricos han sido propuestos por diversos autores (Pousse, Bonilla, Cotton, & Margerin, 2006; Sabetta & Pugliese, 1996). El método propuesto por Pousse, et al. (2006) son modelos de acelerogramas en dominio de tiempo 33

basado en la suposición de que la fase es al azar y que la envolvente de tiempo puede ser representada por una distribución logarítmica normal de las ondas P y S combinadas con una función de algebra-exponencial que representa la envolvente de las ondas P y coda. Además, el contenido de frecuencia de la señal es no estacionaria y sigue un modelo 𝜔-cuadrado modificado. La envolvente de onda propuesto por Pousse, et al. (2006) se observa en la Figura 4.15, esta muestra las diversas partes de la señal: ondas P, ondas S y ondas de coda. Este método depende de cuatro indicadores comunes en la ingeniería sísmica.

Figura 4.15. Envolvente del sismograma sintético propuesto por Pousse, et al. (2006)

4.4.10. El método de Elemento Finito La razón principal por la cual se utiliza elemento finito en el campo de la dinámica es que permite modelar geometrías complejas y leyes constitutivas así como heterogeneidades fuertes (Curnier, 1994; Hughes, 1987), Por otra parte, se adapta bien en la simulación de interacción dinámica suelo-estructura (Clough & Penzien, 1993). 4.4.10.1. Formulación fuerte Para un dominio acotado Ω, la formulación solida se corresponde con el análisis del equilibrio local de una manera exacta, (Salençon, 2001): ∇. 𝜎 + 𝑓 = 𝜌𝑎 𝑒𝑛 Ω

(4.54)

𝜎. 𝑛 = 𝑇 𝑑

(4.55)

𝑒𝑛 𝜕Ω

34

Donde ∇ es el vector gradiente, 𝜎 es el esfuerzo de tensión de Cauchy, 𝑓 el cuerpo de fuerzas, 𝑎 un vector de aceleración y 𝑇 𝑑 fuerzas de superficie. La continuidad de tracciones en cualquier interfaz también debe ser satisfecho: [𝜎] 𝑛 = 𝑡 = 0 𝑒𝑛 ∑

(4.56)

∑.

Donde [𝜎]∑. Denota el salto a través de una superficie de discontinuidades ∑., 𝑛 es el vector normal de ∑. y 𝑡 es el vector de tensión. Teniendo en cuenta el Principio del Trabajo Virtual (Salengon, 2001), el equilibrio ̂ como: global puede ser escrito para cualquier campo de velocidad virtual 𝑈 ̂ )𝑑Ω + ∫ 𝑓. 𝑈 ̂ 𝑑Ω + ∫ 𝑇 𝑑 . 𝑈 ̂ 𝑑a = ∫ 𝜌𝑎. 𝑈 ̂ 𝑑Ω −∫Ω 𝜎: 𝑑̂(𝑈 Ω ∂Ω Ω

(4.57)

Donde 𝑑̂ es la deformación por tensión de velocidad. También puede ser expresado mediante una forma equivalente que implica desplazamientos virtuales en vez de velocidades virtuales (Salençon, 2001): −∫Ω 𝜎: 𝜀̂(𝑢̂)𝑑Ω + ∫Ω 𝑓. 𝑢̂𝑑Ω + ∫∂Ω 𝑇 𝑑 . 𝑢̂𝑑a = ∫Ω 𝜌𝑎. 𝑢̂𝑑Ω

(4.58)

Donde 𝜀̂ es el esfuerzo de tensión. 4.4.10.2. Formulación débil Para un medio isotrópico lineal expuesto a pequeñas deformaciones, la ley constitutiva es escrita como la clásica ley de Hooke: 𝜎 = 𝐶: 𝜀 = 𝜇 (tr𝜀) Π + 2𝜇𝜀

(4.59)

O en términos de desplazamiento: 𝜎𝑖𝑗 = 𝜆(∑3𝑖=1 𝜀𝑖𝑖 )𝛿𝑖𝑗 + 2𝜇𝜀𝑖𝑗

(4.60)

La ley constitutiva puede ser reformulada en la ecuación 4.57 para obtener la forma débil de equilibrio del medio. (Le Tallec, 2007; Salençon, 2001) : 𝜕2 𝑢

−∫Ω 𝜀(𝑢): 𝐶: 𝜀(𝑢̂)𝑑Ω + ∫Ω 𝑓. 𝑢̂𝑑Ω + ∫∂Ω 𝑇 𝑑 . 𝑢̂𝑑a − ∫Ω 𝜌 𝜕𝑡 2 𝑢̂𝑑Ω = 0

35

(4.61)

Para cualquier campo de desplazamiento 𝑢̂ ∈ 𝑆, 𝑆 siendo el espacio de todos los campos cinemáticamente admisibles, y con las siguientes condiciones iniciales en 𝑢: 𝑢(𝑥, 0) = 𝑢0 (𝑥) 𝑦 𝑢̇ (𝑥, 0) = 𝑣0 (𝑥) (𝑥 𝜖 Ω). En lugar de buscar la solución exacta del problema, se puede tratar de estimar una solución aproximada, reduciendo al mínimo el término lateral izquierdo en la ecuación 4.61 denotado 𝑊(𝑢̂) en el espacio de los campos cinematicamente admisibles 𝑢̂ . La ecuación 4.61 se denomina formulación débil del problema. 4.4.10.3. Minimización aproximada: Método de Galerkin La idea principal de la aproximación de Galerkin es reducir al mínimo el término de la izquierda en la ecuación 4.61 no en el espacio de todos los campos cinemáticamente admisibles S pero en un semi-espacio de dimensión finita (d. M. Bonnet & Frangi, 2007; Le Tallec, 2007; Salençon, 2001). La aproximación por tanto consiste en la definición del semi-espacio 𝑆 ℎ , por ejemplo: 𝑆 ℎ = {𝑢̂ℎ = ∑𝑛𝑖=1 𝑎𝑖 𝑁𝑖 + 𝑢𝑑 , 𝑎𝑖 𝜖 ℜ𝑛 } ⊂ 𝑆

(4.62)

Donde 𝑢̂ℎ puede ser expresado a partir de cualquier desplazamiento 𝑢𝑑 satisfaciendo las condiciones limite y n funciones 𝑁𝑖 explícitamente dadas en el espacio de los desplazamientos siendo cero en las fronteras cuando se es prescrito el desplazamiento. Los 𝑛 parámetros escalares 𝑎𝑖 , deberían ser identificados por el siguiente proceso de minimización: min 𝑊(𝑢̂ℎ )

(4.63)

̂ℎ ∈ 𝑆 ℎ 𝑢

Esta formulación se denomina variación semi-discretizada desde que el momento derivado sigue siendo continuo y también tiene que ser aproximado considerando un tiempo discretizado (Hughes, 1987; Reddy, 2006; Zienkiewicz & Taylor, 2005). 4.4.10.4. Discretización del dominio A través del método de Galerkin, el espacio de todos los campos de desplazamiento cinemáticamente admisibles S es aproximado por el semi-espacio 𝑆 ℎ . Una forma práctica de hacerlo consiste en aproximar el dominio exacto Ω por un conjunto de triángulos o áreas cuadriláteras Ω𝑒 (discretizacion espacial) tal que: 𝑛

𝑒 Ω ≃ ∑𝑒=1 Ω𝑒

(4.64)

Estas pequeñas área geométricas se denominan elementos finitos y su dimensión promedio es denotada Δℎ . En el caso de elementos 3D se pueden considerar: tetraedros, pentaedros o hexaedros por ejemplo. (Bathe, 1995; Hughes, 1987). 4.4.11. Varios tipos de elementos finitos 36

Varias órdenes de interpolación. Como se muestra en la figura 4.16, generalmente se considera desplazamientos nodales para definir el campo de desplazamiento en un elemento finito triangular: 𝑢ℎ (𝑥) = 𝑁𝑖 (𝑥)𝑢𝑖 + 𝑁𝑗 (𝑥)𝑢𝑗 + 𝑁𝑘 (𝑥)𝑢𝑘

(4.65)

Las funciones de forma pueden ser lineales, esto es en dos dimensiones: 𝑁1 (𝑥1 , 𝑥2 ) = 𝑎𝑙 + 𝑏𝑙 𝑥1 + 𝑐𝑙 𝑥2

𝑙 = 𝑖, 𝑗, 𝑘

(4.66)

Figura 4.16. Discretización del dominio Ω por elementos finitos Ω𝑒 y una descripción de los desplazamientos nodales (caso 2D). (J.F. Semblat & Pecker, 2009)

También es posible seleccionar funciones con órdenes de segundo grado o mayores (Hughes, 1987; Reddy, 2006; Zienkiewicz & Taylor, 2005). Para un elemento cuadrilátero la forma de la función se observa en la figura 4.17. Las funciones de formas lineales y cuadráticas son definidas como: Elementos lineales cuadriláteros (4 nodos): 1

𝑁1 (𝑟, 𝑠) = 4 (1 − 𝑟)(1 − 𝑠) 1

𝑁2 (𝑟, 𝑠) = 4 (1 + 𝑟)(1 − 𝑠)

(4.67)

1

𝑁3 (𝑟, 𝑠) = 4 (1 + 𝑟)(1 + 𝑠) 1

{𝑁4 (𝑟, 𝑠) = 4 (1 − 𝑟)(1 + 𝑠)

37

Elementos cuadriláteros cuadráticos (8 nodos): 1

𝑁1 (𝑟, 𝑠) = 4 (1 − 𝑟)(1 − 𝑠)(−1 − 𝑟 − 𝑠) 1

𝑁2 (𝑟, 𝑠) = 4 (1 + 𝑟)(1 − 𝑠)(−1 + 𝑟 − 𝑠) 1

𝑁3 (𝑟, 𝑠) = 4 (1 + 𝑟)(1 + 𝑠)(−1 + 𝑟 + 𝑠) 1

𝑁4 (𝑟, 𝑠) = 4 (1 − 𝑟)(1 + 𝑠)(−1 − 𝑟 + 𝑠)

(4.68)

1

𝑁5 (𝑟, 𝑠) = 4 (1 − 𝑟 2 )(1 + 𝑠) 1

𝑁6 (𝑟, 𝑠) = 4 (1 − 𝑠 2 )(1 + 𝑟) 1

𝑁7 (𝑟, 𝑠) = (1 − 𝑟 2 )(1 + 𝑠) 4 1

2 { 𝑁8 (𝑟, 𝑠) = 4 (1 − 𝑠 )(1 − 𝑟)

Como en la figura 4.17, cada función de forma es cero en todos los nodos. Debido al mayor grado de las funciones en el caso cuadrático (parte inferior), es posible recuperar las variaciones más rápido en un elemento cuadrático que un elemento lineal. Esta observación es particularmente importante para la simulación de propagación de onda.

Figura 4.17. Funciones de forma lineal (arriba) y cuadrática (abajo) elementos finitos. (J.F. Semblat & Pecker, 2009)

Elementos finitos de referencia: Además de la interpolación del campo de desplazamiento dentro de cada uno de los elementos finitos es conveniente considerar un elemento de referencia, fijando geometría normalizada, para llevar a 38

cabo cálculos básicos siempre de la misma manera. (Bathe, 1996; Hughes, 1987), Como se muestra en la figura, el elemento cuadrilátero de referencia es un cuadrado y cualquier geometría cuadrilátera puede ser definida por una simple transformación entre elementos reales y el de referencia. La matriz Jacobiana de esta transformación debe ser calculada desde los elementos integrales de la misma manera. Por cada uno de los elementos finitos, por lo que se hace necesario interpolar ambos el campo de desplazamiento y la geometría. Una manera más simple de realizar ambas interpolaciones es el de considerar las mismas funciones de forma; tales elementos finitos se denominan elementos isoparamétricos (Reddy, 2006; Zienkiewicz & Taylor, 2005) (Figura 4.18).

Figura 4.18. Coordenadas intrínsecas de los elementos finitos cuadriláteros distorsionados: elementos distorsionados (arriba) se definen a partir de elementos de referencia (abajo). (J.F. Semblat & Pecker, 2009).

39

Figura 4.19. Varios tipos de elementos finitos de diferente orden p. (Hughes, 1987)

Los diferentes tipos de elementos isoparamétricos triangulares y cuadriláteros, en varios ordenes p (Hughes, 1987) son mostrados en la figura 4.19. Lineal (arriba) y cuadrática (centro), también se muestran casos de orden superior, un elemento triangular de orden 3 (parte inferior izquierda) y un elemento cuadrilátero de 4 orden (parte inferior derecha). 4.4.12. Algoritmos de integración del tiempo La formulación variacional semi-discreta solo se refiere a la discretización espacial. Es también necesario considerar una discretización del tiempo para aproximar las derivadas de tiempo. Por lo tanto se considera 𝑚 + 1 tiempos discretos 𝑡𝑛 tales 𝑡𝑚𝑎𝑥 como 𝑡𝑛 = 𝑛Δ𝑡, 𝑛 = 0, 𝑚. El paso de tiempo ∆𝑡 es definido como ∆𝑡 = 𝑚 donde 𝑡𝑚𝑎𝑥 es el máximo tiempo de la simulación. Comenzando a partir de condiciones iniciales, {{𝑈0 }, {𝑈̇0 }} El algoritmo de integración de tiempo procede a aproximar soluciones en cada paso de tiempo nombrados ({𝑈𝑛 }, {𝑈̇𝑛 }, {𝑈̈𝑛 }). Varios esquemas de integración de tiempo están disponibles: explícitos, implícitos, el algoritmo de múltiples etapas, etc (Hughes, 1987; Zienkiewicz & Taylor, 2005). A continuación solo se discutirá algunos esquemas de integración de tiempo clásico y sus propiedades (ventajas e inconvenientes) en el campo de la dinámica y propagación de onda.

40

4.4.12.1. El método de diferencia central El método de diferencia central es un esquema de integración temporal explicita simple (Bathe, 1996; Reddy, 2006). La velocidad y la aceleración son expresas de forma sencilla: 1 {𝑈̇𝑛+1 } = 2∆𝑡 ({𝑈𝑛+1 } − {𝑈𝑛−1 })

(4.69)

1 {𝑈̈𝑛 } = ∆𝑡 2 ({𝑈𝑛+1 } − 2{𝑈𝑛 } + {𝑈𝑛−1 })

(4.70)

Por tanto, se produce una aproximación de segunda orden clásica del vector de aceleración, Si se tiene en cuenta una matriz de masa concentrada, la expresión del algoritmo de diferencia central es totalmente explícita. Sin embargo este algoritmo es solo condicionalmente estable (J.F. Semblat & Pecker, 2009). 4.4.12.2. El método de Newmark El método de Newmark es un enfoque preciso de segundo orden ampliamente utilizado en el campo de la dinámica (Hughes, 1987; Zienkiewicz & Taylor, 2005). Se trata de un algoritmo de dos parámetros (𝛽, 𝛾) definidos de la siguiente manera: [𝑀]{𝑈̈𝑛+1 } + [𝐶]{𝑈̇𝑛+1 } + [𝐾]{𝑈𝑛+1 } = {𝐹𝑛+1 } ∆𝑡 2 {𝑈𝑛+1 } = {𝑈𝑛 } + ∆𝑡{𝑈̇𝑛 } + [(1 − 2𝛽){𝑈̈𝑛 } + 2𝛽{𝑈̈𝑛+1 }] 2

(4.71)

{𝑈̇𝑛+1 } = {𝑈̇𝑛 } + Δ𝑡[(1 − 𝛾){𝑈̈𝑛 } + 𝛾{𝑈̈𝑛+1 }] Las condiciones de estabilidad del algoritmo son las siguientes: Estabilidad incondicional 1 2

≤ 𝛾 ≤ 2𝛽

(4.72)

Estabilidad condicional 1

𝛾

𝛾 ≥ 2 ; 𝛽 < 2 ; 𝑤 ℎ ∆𝑡 ≤ Ω𝑐𝑟𝑖𝑡

(4.73)

Donde 𝑤 ℎ es la frecuencia circular aproximada de elementos finitos Δℎ 𝑦 Ω𝑐𝑟𝑖𝑡 la frecuencia critica:

Ω𝑐𝑟𝑖𝑡 =

1 𝛾 1 𝜉(𝛾− )+√ −𝛽+𝜉 2 (𝛾− ) 2

2 𝛾 −𝛽 2

2

2

(4.74)

Donde 𝜉 es el coeficiente de amortiguamiento 41

Opciones clásicas para los parámetros del algoritmo son:   

1

1

Método de aceleración promedio: Los parámetros son 𝛽 = 4 𝑦 𝛾 = 2. El algoritmo es por tanto implícito e incondicionalmente estable. 1 1 Método de aceleración lineal: Los parámetros son 𝛽 = 6 𝑦 𝛾 = 2. El algoritmo es implícito y condicionalmente estable con Ω𝑐𝑟𝑖𝑡 = 2√3 1 Método de diferencia central: Seleccionando 𝛽 = 0 𝑦 𝛾 = 2 en el esquema de Newmark conduce al método de diferencia central. El algoritmo es explícito y condicionalmente estable con Ω𝑐𝑟𝑖𝑡 = 2.

Para otros valores de los parámetros 𝛽 𝑦 𝛾, las condiciones de estabilidad son mostradas en la figura con respecto a 𝐴1 𝑦 𝐴2 definidas como (Curnier, 1994; Hughes, 1987): 𝐴1 = 1 − 𝐴2 = 1 −

Ω2 1 (𝛾+ )] 2 2 1+2𝛾𝜉Ω+𝛽Ω2

[𝜉Ω+

(4.75)

1 2 1+2𝛾𝜉Ω+𝛽Ω2

[2𝜉Ω+Ω2 (𝛾− )]

(4.76)

Con: Ω = 𝑤 ℎ Δ𝑡 1

Fuera del triangulo (𝛾 < 2 𝑦 2𝛽 − 𝛾 < 0), En el esquema de Newmark es inestable 1

(Curnier, 1994). En el triángulo puede ser incondicionalmente estable (𝛾 ≥ 2 𝑦 2𝛽 − 𝛾 < 0). Los rangos correspondientes para el radio espectral 𝜌 (Curnier, 1994; Reddy, 2006) se dan en la figura 4.20. 4.4.12.3. Método de pasos múltiples Entre los métodos de pasos múltiples esta, El método 𝛼 propuesto por Hughes (1987) el cual tiene la siguiente forma: [𝑀]{𝑈̈𝑛+1 } + (1 + 𝛼)[𝐶]{𝑈̇𝑛+1 } − 𝛼[𝐶]{𝑈̇𝑛 } + (1 + 𝛼)[𝐾]{𝑈𝑛+1 } − 𝛼[𝐾]{𝑈𝑛 } = {𝐹(𝑡𝑛+𝑎) } 𝑐𝑜𝑛 𝑡𝑛+𝑎 = (1 + 𝛼)𝑡𝑛+1 − 𝛼𝑡𝑛 = 𝑡𝑛+1 + 𝛼Δ𝑡 (4.77)

42

Figura 4.20. Condiciones de estabilidad de integración-tiempo de Newmark (Curnier, 1994; Hughes, 1987).

El cual corresponde al algoritmo de Newmark para 𝛼 = 0 Este método es a menudo llamado el 𝑚𝑒𝑡𝑜𝑑𝑜 𝛼 − 𝐻𝐻𝑇. Esta incondicionalmente 1

(1−𝛼)2

1−2𝛼

estable cuando − 3 < 𝛼 < 0, 𝛾 = 2 𝑦 𝛽 = 4 . Su interés principal es que puede suprimir fácilmente componentes de alta frecuencia (Usando algoritmos de amortiguamiento). Otro método de este tipo fue propuesto por Wood et al. (1981). En este método únicamente el término inicial es modificado: (1 − 𝛼𝐵 )[𝑀]{𝑈̈𝑛+1 } + 𝛼𝐵 [𝑀]{𝑈̈𝑛 } + [𝐶]{𝑈̇𝑛+1 } + [𝐾]{𝑈𝑛+1 } = {𝐹𝑛+1 } 1

𝛾

(4.78) 1

1

Este método es estable incondicionalmente para 𝛼𝐵 ≤ 2 , 𝛽 ≥ 2 ≥ 4 , 𝛼𝐵 + 𝛾 ≥ 2. También es posible considerar esquemas de tiempo de integración de orden superior, utilizando por ejemplo procedimientos paso a paso. La exactitud de las respuestas dinámicas calculadas parecen ser mejoradas al usar tales procedimientos (Fung, 1997). Varios de otros esquemas de integración han sido propuestos en el esquema de la dinámica: Wilson, Park, Houbolt, etc. (Bathe, 1996; Hughes, 1987; Zienkiewicz & Taylor, 2005).

4.4.13. Modelado de propagación de ondas en medios infinitos El Modelado de Propagación de onda en medios semi-infinitos puede ser difícil dependiendo del método numérico considerado. Los métodos integrales de la 43

ecuación de contorno permiten una descripción precisa de la propagación de ondas en medios infinito (M. Bonnet, 1999). Los métodos de elementos finitos o espectrales, son necesario para evitar reflexiones de la onda artificiales en los límites de malla. Varias técnicas que resuelven este problema se discutirá a continuación: contornos absorbentes, elementos infinitos y capas absorbentes (Festa & Nielsen, 2003; Givoli, 1992; Meza-Fajardo & Papageorgiou, 2008; Modaressi & Benzenati, 1992). 4.4.13.1. Contornos Absorbentes en 2D En los casos de 2D, la solución teórica para las condiciones óptimas de absorción no es tan simple como en el caso 1D. Se considerará un ejemplo en 2D para el diseño de absorción de las condiciones de contorno: un dominio cilíndrico es modelado por un cuarto de un disco discretizado con elementos finitos (deformación plana). Como se muestra en la Figura 4.21, elementos discretos están conectados en el límite del dominio. El objetivo es tener una determinación teórica de las características óptimas de los elementos de absorción con el fin de minimizar o eliminar por completo reflexiones de la onda espurias. Este ejemplo fue propuesto por Dubreucq y Piau para analizar la interacción sueloestructura. La primera etapa consiste en la obtención de la ecuación de propagación para las ondas cilíndricas (Bisch, Langeoire, Prat, & Semblat, 1999). El medio se supone elástico, lineal e isótropo en el marco de pequeñas deformaciones. El campo de desplazamiento se supone que depende sólo de la distancia radial r.

Figura 4.21. Malla en 2D con amortiguador y resortes en su contorno (Bisch, et al., 1999)

44

𝜆

𝜇

Los términos 𝑅 y − 𝑅 corresponden a las rigideces de línea para cada tipo de onda. 𝑒

𝑒

También es necesario para determinar los rasgos característicos de los resortes y amortiguadores en el marco de una formulación de elementos finitos (dependiendo del tipo de elemento) y tener en cuenta la discretización de la frontera (Bisch, et al., 1999). Para una discretización dada de elementos finitos y de excitación (ondas P o S), es por lo tanto posible determinar la rigidez óptima y los valores de atenuación para eliminar las reflexiones falsas en los límites del dominio. No obstante, se debe notar que estos valores óptimos se determinan suponiendo que el radio del dominio cilíndrico 𝑅𝑒 , debe ser grande y por lo tanto obtiene la simplificación de los componentes de la tensión.

Figura 4.22. Esquema de elementos infinitos en 2D a lo largo de una dirección (arriba) o de dos direcciones (abajo), propuesto por J.F. Semblat and Pecker (2009).

4.4.13.2. Caso bidimensional En 2D, el principio del método es el mismo que en 1D y en 2D el elemento puede ser infinito a lo largo de una o dos direcciones. Figura 4.22 (parte superior) muestra un elemento infinito 2D a lo largo de una sola dirección. En este caso, Chadwick, Bettess, and Laghrouche (1999) propusieron un elemento radialmente infinito con las siguientes funciones de mapeo: 2𝑠

1

𝑀1 = (− 1−𝑠) (− 2 𝑟(1 − 𝑟)) 1+𝑠

1

𝑀2 = (1−𝑠) (− 2 𝑟(1 − 𝑟)) 2𝑠

𝑀3 = (− 1−𝑠) (1 + 𝑟)(1 − 𝑟)

(4.79) 45

1+𝑠

𝑀4 = (1−𝑠) (1 + 𝑟)(1 − 𝑟) 2𝑠

1

𝑀5 = (− 1−𝑠) 2 𝑟(1 − 𝑟) 1+𝑠 1

𝑀2 = (1−𝑠) 2 𝑟(1 − 𝑟) Los elementos infinitos representados en la Figura 4.22 (arriba) pueden estar relacionados con elementos finitos estándar a través de estas transformaciones. En el caso de un elemento infinito a lo largo de dos direcciones (Figura 4.22, parte inferior), las expresiones anteriores pueden ser fácilmente generalizadas (Bettess & Dasgupta, 1992). 4.4.14. Efecto de sitio: análisis 1D, 2D Como es meritorio saber las amplificaciones de las ondas sísmicas están influenciadas por las propiedades mecánicas de los estratos. De hecho, el cambio de velocidad entre las capas del suelo (heterogeneidades verticales) modifica la amplificación de movimiento en la superficie libre. Como se representa en la Figura 4.23 (a la izquierda), que es principalmente un efecto 1D. La geometría del depósito es también un factor determinante (Figura 4.23, derecha). Se puede estar caracterizado por una profundidad media o de un formulario detallado para una cuenca sedimentaria. En este último caso (heterogeneidades laterales), las ondas sísmicas son atrapados en la cuenca que conduce a una amplificación mayor que a través de la suposición de 1D. En el caso en el que sólo se considera la profundidad media de las capas (caso unidimensional), el factor de amplificación se puede estimar teóricamente. La estimación del movimiento del suelo en los códigos se lleva a cabo principalmente a través de esos métodos simplificados (Pitilakis, Raptakis, & Makra, 1999). Considerando que, cuando las heterogeneidades laterales son grandes, es necesario investigar la propagación de ondas sísmicas en dos o tres dimensiones (Bard & Bouchon, 1985).

Figura. 4.23. Amplificación de las ondas sísmicas: fenómenos principales en 1D (izquierda) y estructuras geológicas 2D/3D (derecha). (J.F. Semblat & Pecker, 2009)

4.4.15. Región de absorción Las regiones de absorción de energía se utilizan para evitar los reflejos en los límites de simulaciones de propagación de ondas con elementos finitos. Para los modelos 46

de dominio de frecuencia es importante utilizar excelentes contornos absorbentes para simular un espacio infinito evitando posibles artefactos de reflexiones en los límites. Según Drozdz (2008) asegura que en el proceso de modelado con software comerciales basados en elementos finitos es imprescindible la aplicación de regiones absorbentes viscoelásticas y capas perfectamente emparejadas, ya que permiten una reducción en la geometría del modelo. En este enfoque, el factor de pérdidas (η) del módulo complejo viscoelástico se aumenta gradualmente de 0 a 1 en la región de absorción. En comparación con capas perfectamente emparejadas, este enfoque requiere más elementos que son más robusto y fácil de usar ya que la absorción de más de una distancia fija no depende de la dirección de la propagación de frentes de onda, es decir, que también absorben las ondas que tienen velocidades de energía y la fase de signos opuestos (Ryden & Castaings, 2009). En un caso simétrico axial 2D, la siguiente ecuación se utiliza para definir el módulo complejo en r-dirección dentro de la región de absorción y se da como: 𝑟−𝑟𝑎 3

𝐸 = 𝐸𝑟𝑒𝑎𝑙 [1 + 𝑖 (𝑅

𝑎𝑏𝑠𝑟

) ]

(4.80)

Donde 𝑟 es la posición a lo largo de la dirección y 𝑟𝑎 y 𝑅𝑎𝑏𝑠𝑟 son la posición inicial y la longitud respectivamente. 4.4.16. Geometría adaptable y tamaño de la malla Ryden and Castaings (2009) proponen desarrollar un modelo de elementos finitos para la simulación eficiente de onda de superficie utilizando una geometría adaptativa y un tamaño de malla para cada frecuencia. El espacio de modelado se divide en una región analizada y una región de absorción. La región analizada es el espacio de interés y la región de absorción sólo se utiliza para simular condiciones de borde infinito. La región analizada se mantiene fija salvo el tamaño de cada elemento y las regiones absorbentes están optimizados para cada frecuencia modelada. Esto significa que el modelo total y el tamaño de los elementos son grande cuando se resuelve para las frecuencias bajas y pequeñas, que para cuando se resuelve para las frecuencias altas.

47

Figura 4.24. Esquema del modelo de elementos finitos adaptativo con regiones analizadas y de absorción aplicando el análisis de estructuras de pavimentos. (Ryden & Castaings, 2009)

Ryden and Castaings (2009) citado por Obando, Park, Ryden, and Ulriksen (2010) sugieren que el modelado debe realizarse utilizando un enfoque de elementos finitos adaptativo en el dominio de frecuencia (0-10 Hz) usando regiones absorbentes viscoelásticas para simular límites infinitos; además que la región de absorción debe extenderse a ambos lados del modelo hasta 3 veces la longitud de onda máxima de cada capa en la dirección horizontal y 5 veces la longitud de onda máxima del semiespacio en la dirección vertical (Figura 4.24). Aplicando la ecuación 4.48, se calcula la longitud de onda máxima (𝜆): 𝜆=𝑓

𝑉𝑝

(4.81)

𝑚𝑖𝑛

Teniendo 𝜆 se diseña la región de absorción 𝐿′ y 𝐻 ′ a como se describe a continuación: 𝐿′ = 3 ∗ 𝜆

(4.82)

𝐻′ = 5 ∗ 𝜆

(4.83)

4.5. ESPECTRO DE RESPUESTA En ingeniería sísmica, el espectro de respuesta brinda un significado conveniente a la sumatoria de respuestas pico de todos los posibles sistemas simples sujeto a un componente particular del movimiento del suelo. Una gráfica de valores pico de respuesta de una cantidad como función del periodo natural de vibración del sistema o cualquier parámetro relacionado como frecuencia 48

angular (𝜔𝑛 ) o Frecuencia natural (𝑓𝑛 ) es denominado espectro de respuesta para esa cantidad. Existen diversos tipos de espectros de respuesta en dependencia de la reacción que se desea comparar, espectros de respuesta de deformación, espectros de respuesta de Pseudo-Velocidad y el más habitual en cálculos sísmicos es el espectro elástico de respuesta, este último relaciona la aceleración, y es de considerable importancia ya que aplicando la segunda ley de Newton es posible conocer las fuerzas estáticas equivalente que soporta una estructura mediante el producto de la aceleración y la masa que soporta la estructura. (Chopra, 2007) 4.5.1. Espectro de Respuesta de Deformación. Este espectro consiste en una gráfica la cual posee en el eje de las ordenadas la deformación (𝑢0 ), en el eje de las abscisas el Periodo natural (𝑇𝑛 ) para un amortiguamiento (𝜉) fijo. La figura 4.25 permite apreciar el procedimiento para la obtención del espectro, dicho espectro es desarrollado para el movimiento sísmico de El Centro, Figura 4.25 (a). La variación de la deformación inducida por el movimiento del suelo es mostrada en la figura 4.25 (b). Para cada sistema el valor pico de deformación es determinado del histograma de deformación. El valor de la amplitud 𝑢0 determinado para cada sistema provee una coordenada en el espectro de respuesta de deformación. Repitiendo estos cálculos para un rango de valores de 𝑇𝑛 mientras 𝜉 permanece constante, se provee el espectro de deformación. Figura 4.25 (c)

Figura 4.25. (a) Aceleración del suelo, (b) Respuesta de deformación de tres sistemas simples con 𝜉 = 2% 𝑦 𝑇𝑛 = 0.5; 1; 2 𝑠𝑒𝑔, (c) Espectro de Respuesta de Deformación para 𝜉 = 2% (Chopra, 2007)

49

4.5.2. Espectro de Respuesta de Pseudo-Velocidad Considera una cantidad V para un sistema simple con una frecuencia natural (𝑤𝑛 ) relacionado con su deformación pico 𝐷 = 𝑢0 debido al movimiento del suelo por: 𝑉 = 𝜔𝑛 ∗ 𝐷 =

2𝜋 𝑇𝑛

∗𝐷

(4.84)

Donde V es llamada Pseudo-velocidad pico. Un ejemplo de este grafico se muestra en la figura 4.26.

Figura 4.26. Espectro de Respuesta de Pseudo-Velocidad. (Chopra, 2007)

4.5.3. Espectro de Respuesta de Pseudo- Aceleración Considera una cantidad A para un sistema simple con una frecuencia natural (𝑤𝑛 ) relacionado con su deformación pico 𝐷 = 𝑢0 debido al movimiento del suelo por: 2𝜋 2

𝐴 = (𝜔𝑛 )2 ∗ 𝐷 = ( 𝑇 ) ∗ 𝐷

(4.85)

𝑛

Donde A es llamada Pseudo-aceleración pico. Un ejemplo de este grafico se muestra en la figura 4.27.

50

Figura 4.27. Espectro de Respuesta de Pseudo- Aceleración. (Chopra, 2007)

4.5.4. Espectro de Respuesta Combinado D-V-A

Figura 4.28. Espectro de Respuesta Combinado D-V-A para el sismo de El Centro, 𝜉 = 2% (Chopra, 2007)

Los tres espectros proveen directamente cantidades físicas significativas, es por esta razón su necesidad, El espectro de deformación brinda la deformación pico del sistema; el espectro de Pseudo-Velocidad se relaciona directamente con la energía 51

pico almacenada en el sistema durante un sismo, el espectro de PseudoAceleración se relaciona directamente con el valor pico de la fuerza estática equivalente y el cortante basal. Para propósitos prácticos de diseño las tres cantidades espectrales pueden ser representadas en un solo grafico en un papel tetra logarítmico. Todo esto gracias a la siguiente interrelación (Figura 4.28). 𝑇𝑛 2𝜋

∗𝐴=𝑉 =

2𝜋 𝑇𝑛

∗𝐷

(4.86)

52

5. PROCESAMIENTO Y ANÁLISIS DE DATOS

En este apartado se describen los procedimientos para el procesamiento y análisis de los datos, así como también las herramientas seleccionadas para su ejecución. Primeramente se plantean los perfiles 2D basado en estudios previos (Castillo & Zepeda, 2013), y los criterios para la selección de las zonas de transición. Luego se presenta el procedimiento para el modelamiento en elementos finitos codificado en el programa Comsol Multiphysics. Finalmente se explican el esquema de posprocesamiento para el análisis de las funciones de transferencia sísmica así como el cálculo de acelerogramas de superficie y el posterior cálculo de espectros elásticos de respuesta. 5.1. DETERMINACIÓN DE MODELOS DE VELOCIDAD Los registros utilizados para plantear perfiles que contengan las zonas de transición de una zona a otra fueron obtenidos de la tesis monográfica de los Ingenieros Castillo and Zepeda (2013) quienes propusieron 5 zonas para la ciudad de Managua en dependencia de los registros de frecuencia obtenidos de los análisis de razones espectrales H/V. Estos modelos fueron derivados usando de referencia los obtenidos por Escorcia and Ochoa (2013) los cuales se basaron en resultados previamente obtenidos por Faccioli, et al. (1973) y Obando Hernández (2011). La ubicación de las zonas de transición localizadas al noroeste del área urbana de Managua ilustradas como puntos (figura 5.1) corresponden a los límites de zona, que representan una longitud determinada para los modelos de velocidad.

53

Figura 5.1. Mapa de Ubicación de las zonas de transición seleccionadas

54

Definido los puntos que limitan cada una de las zonas en estudio, se proponen los modelos de velocidad de cortante procesados en el software Microsoft Excel 2013 (figura 5.2), que muestran significativamente la variación en profundidad a una velocidad constante que circula por nueve zonas de transición diferentes.

Figura 5.2. Modelo de Velocidad de cortante

5.2. MODELAMIENTO EN DOS DIMENSIONES Para la evaluación de respuesta de modelos de suelo en 2 dimensiones se hará uso del software basado en la teoría de elementos finitos codificado en Comsol Multiphysics 4.2a. Para el análisis de propagación de ondas se usa el módulo de mecánica estructural que permite modelar ondas de corte en dominio de frecuencia asumiendo un medio lineal elástico (figura 5.3). Con este software se podrá modelar el efecto de las zonas de transición en la estimación de la respuesta de sitio representada en términos de espectros elásticos de respuesta.

55

Figura 5.3. Interfaz gráfica de software Comsol Multiphysics para el procesamiento en mecánica de sólidos.

Una vez definidos los modelos de velocidades en 2D, el modelamiento en Comsol Multiphysics se continúa de la siguiente manera:     

Inicialmente se seleccionara el módulo de mecánica estructural para modelamiento de esfuerzos; para luego definir la geometría del modelo que incluye las zonas de transición. Posteriormente se definen geometrías que representan dos estratos finitos de suelo, el primero horizontalmente infinito, mientras el segundo representa la zona de absorción o el espacio semi-infinito en profundidad del suelo. Seguidamente a cada estrato se le asignan propiedades que definen su entrada al modelo lineal elástico (Velocidad, módulo de Young, relación de Poisson, densidad, etc.) Se asigna una fuente excitadora la que se representa como una carga unitaria horizontal para simular la propagación vertical de ondas SH incidente. Finalmente se modela la respuesta del modelo de suelo en dominio de frecuencia en un rango de 1 a 10 Hz. Los resultados son exportados a Matlab para procesamientos adicionales. (Figura 5.4)

56

5.3. CÁLCULO DE RESPUESTA SÍSMICA DE SITIO Una vez completado el modelamiento en Comsol se exportaran los resultados a formato “.mat” de Matlab para ser pos-procesados y así realizar el cálculo de los cocientes espectrales a partir de los espectros de amplitud.

Figura 5.4. Post-procesamiento en rutinas de Matlab describiendo una función de trasferencia

Los resultados obtenidos en Comsol serán luego validados, en primera instancia con las respuestas lineales elásticas en una dimensión obtenida a partir de la teoría de multi-reflección codificada en Matlab (Figura 5.4) (Obando, et al., 2010).

57

(a)

(b)

Figura 5.5. (a) Rutina en Matlab para cálculo de función de transferencia en condiciones de horizontalidad del suelo (b) función de transferencia generada por algoritmo.

La comparación de los modelamientos y el análisis en 1D (Figura 5.5) se hace inicialmente en término de función de transferencia sísmica para el mismo rango de frecuencias. 5.3.1. Cálculo de acelerogramas sintéticos en superficie La respuesta sísmica de sitio se lleva a cabo calculando la respuesta del suelo en término de espectros elásticos de respuesta. Para ello se usara el registro deconvolucionado del evento principal del terremoto de Managua de 1972. 5.3.2. Proceso de convolución Para el cálculo de un registro en superficie debido al registro de entrada se hará uso del proceso de convolución en dominio de frecuencia. Por lo tanto la convolución se hace tomando la función de transferencia compleja modelada en Comsol y la transformada de Fourier del registro de entrada. Luego mediante la inversa de la transformada de Fourier se calcula la respuesta teórica del modelo en superficie libre, este procedimiento se hace un una rutina en Matlab diseñada específicamente para este propósito (Figura 5.6).

58

Figura 5.6. (a) Algoritmo en Matlab para generar el acelerograma en superficie (b) función de transferencia sísmica y acelerograma en basamento y en superficie.

5.3.3. Cálculo de espectros de respuesta Una vez obtenidos los acelerogramas en superficie se calculan los espectros de respuesta de Deformación, Pseudo-velocidad y Pseudo-aceleración mediante un algoritmo basado en la teoría del método de diferencia central en dominio de tiempo descrita por Chopra (2007). Los espectros elásticos calculados para cada modelo se obtendrán antes, durante y después de la transición y así analizar la variabilidad de las aceleraciones espectrales en zonas con variación en la horizontalidad de suelo marcadas.

59

Figura 5.7. (Izquierda) Post-procesamiento en Matlab para generar espectros de respuesta (derecha superior) espectro de deformación, (derecha central) espectro de velocidad y (derecha inferior) espectro de aceleración

El cálculo de los espectros se hace en Matlab usando el algoritmo de diferencia central en dominio de tiempo. En la figura 5.7 se muestra las salidas del cálculo el cual ilustra los espectros de desplazamiento, Pseudo-velocidad y Pseudoaceleración.

60

6. RESULTADOS

En este capítulo se muestran los resultados generados según el esquema de procesamiento expuesto en el apartado anterior. La primera parte corresponde al planteamiento de los 9 modelos en 2D con variación lateral en las zonas propuestas por Castillo and Zepeda (2013) en la ciudad de Managua, haciendo el trazado de líneas en la transición de las zonas para ser evaluados posteriormente. Luego se muestran las funciones de transferencia sísmica generadas por los modelos considerando 3 condiciones, la primera y la tercera establecen que los estratos son horizontalmente infinitos llamándoles para cada caso antes de la transición (A.T) y después de la transición (D.T), y la segunda condición considerando la heterogeneidad del suelo lo define como el lugar de la transición (T), los gráficos son superpuestos para apreciar las posibles variaciones en contenido de frecuencia y amplitud. Por otra parte se determinaron los acelerogramas en superficie obtenidos a partir de la convolución entre las funciones de transferencia sísmica y registro convertido a roca del terremoto de 1972 obtenida mediante deconvolucion por Castillo and Zepeda (2013). Finalmente se generan los espectros elásticos de respuesta de deformación, Pseudo-velocidad y Pseudo-aceleración superpuestos para las tres condiciones antes mencionadas.

6.1. MODELOS EN 2D La creación de los modelos en 2D es el principio para el análisis del efecto de la no horizontalidad del suelo. Se elaboraron a partir de la propuesta de zonificación presentada por Castillo and Zepeda (2013), quienes plantearon 5 zonas en las cuales definen modelos de velocidad de cortante promedio con 3 estratos finitos. En la figura 6.1 se muestra la selección de 9 zonas de transición donde se dan los cambios de suelo según el mapa propuesto por Castillo and Zepeda (2013). Los modelos en 2D propuestos se presentan en un plano cartesiano, los datos de las abscisas se obtienen de la vista en planta proporcionada por el mapa, y el de las ordenadas se obtiene de las columnas estratigráficas.

61

Figura 6.1. Mapa propuesto por Castillo and Zepeda (2013) mostrando la distribución de zonas. Líneas de transición para creación de perfiles 2D generada por los autores.

62

Las velocidades de corte de cada estrato finito se promediaron obteniéndose de esa manera un Vs representativo a lo largo del perfil y para cada estrato (Tabla 6.1).

Figura 6.2. Perfiles Estratigráficos (2D): (a) Bo. El Bóer, (b) Bo. Jorge Dimitrov, (c) Bo. Recreo Sur.

Se muestra en el Modelo 1 (Figura 6.2a) los máximos contrastes de impedancia respecto a la profundidad de la roca elástica alcanzando una depresión de hasta 38 m. En el caso de la Figura 6.2b (Modelo 2) la variación lateral del suelo es más suavizada pasando de una profundidad de 16 a 10 m aproximadamente. En otra instancia se ilustra el Modelo 6 (Figura 6.2c), siendo uno de los más lineales

63

obtenidos luego del post-procesamiento, donde transita a profundidades mayores a 10 m, perteneciendo a frecuencias bajas menores a 6 Hz. De manera general se puede decir que 5 perfiles presentan pendientes pronunciadas mayores al 10% obtenidas por la relación entre estratos más profundos a más someros, 4 perfiles son más suavizados con pendientes menores al 10%, y un perfil ligeramente mayor al 10% (Tabla 6.2). El resto de los modelos pertenecientes a los perfiles 2D se muestran en el Anexo 2. La Tabla 6.1 agrupa todos los modelos de perfiles 2D generados, se definen 3 estratos finitos y un semi-espacio extraídos de los modelos de velocidad de cortante propuestos por Castillo and Zepeda (2013); por otra parte como cada modelo abarca entre 3 a 2 zonas representativos a las líneas de transición, se hace un promedio de las velocidades y el resultado es el que se muestra en la tabla.

64

Tabla 6.1. Agrupación de modelos y velocidades promedios para cada estrato

Modelos Modelo 1 (Bo. El Bóer)

Estrato H (m) Vs (m/s) 1 Zona 1 12.50 Zona 2 3.81 Zona 4 3.72 187.61 2 7.77 7.23 2.92 330.15 3 18.10 9.24 3.33 374.21 Semi-espacio 797.67 Modelo 2 1 Zona 2 3.81 Zona 4 3.72 Zona 3 4.33 188.09 (Bo. William 2 7.23 2.92 7.34 336.20 Díaz) 3 9.24 3.33 5.00 352.02 Semi-espacio 857.78 Modelo 3 1 Zona 3 4.33 Zona 4 3.72 Zona 5 3.95 194.24 (Bo. Jorge 2 7.34 2.92 3.50 380.27 Dimitrov) 3 5.00 3.33 3.25 387.39 Semi-espacio 929.05 Modelo 4 1 Zona 1 12.50 Zona 2 3.81 197.83 (Bo. El Cortijo) 2 7.77 7.23 329.23 3 18.10 9.24 388.82 Semi-espacio 816.50 Modelo 5 1 Zona 1 12.50 Zona 4 3.72 193.76 ( M.H.C.P) 2 7.77 2.92 339.39 3 18.10 3.33 399.67 Semi-espacio 816.50 Modelo 6 1 Zona 2 3.81 Zona 3 4.33 198.55 (Bo. Recreo Sur) 2 7.23 7.34 338.30 3 9.24 5.00 355.54 Semi-espacio 906.67 Modelo 7 1 Zona 4 3.72 Zona 5 3.95 180.47 (Bo. San José 2 2.92 3.50 387.96 Oriental) 3 3.33 3.25 387.20 Semi-espacio 866.91 Modelo 8 1 Zona 3 4.33 Zona 5 3.95 207.77 (Bo. La 2 7.34 3.50 404.41 Esperanza) 3 5.00 3.25 408.59 Semi-espacio 1013.57 Modelo 9 1 Zona 5 3.95 Zona 1 12.50 207.05 (Planes de 2 3.50 7.77 395.35 Altamira) 3 3.25 18.10 441.87 Semi-espacio 923.41 65

Es importante recalcar que los 3 primeros modelos presentan 3 zonas que limitan la línea de transición, en este caso la longitud del modelo estriba alrededor de los 120 m, mientras el resto de modelos que limitan solo dos zonas la longitud se aproxima a los 60 m. Tabla 6.2. Resumen de pendientes (%) para cada zona de transición

MODELOS

PROFUNDIDAD MAXIMA (m)

PROFUNDIDAD MINIMA (m)

MODELO 1 MODELO 2 MODELO 3 MODELO 4 MODELO 5 MODELO 6 MODELO 7 MODELO 8 MODELO 9

38.37 20.28 16.67 38.37 38.37 20.28 10.70 16.67 38.37

9.97 9.97 9.97 20.28 9.97 16.67 9.97 10.70 10.70

DISTANCIA ENTRE P.MAX Y P.MIN (m) 113.14 61.07 61.07 61.07 61.07 61.07 61.07 61.07 61.07

PENDIENTE (%)

25.10 16.88 10.97 29.62 46.50 5.91 1.20 9.77 45.30

6.2. FUNCIÓN DE TRANSFERENCIA SÍSMICA 6.2.1. Validación del uso de elementos finitos para propagación de ondas 1D Luego de proponer los modelos de perfiles estratigráficos, es requerido determinar si las funciones de transferencia sísmica obtenidas mediante el procesamiento en Comsol y post-procesadas en Matlab son aproximadamente iguales en rango frecuencial, a las obtenidas mediante el algoritmo en Matlab basado en la teoría de propagación de ondas unidimensionales elaborado por Robinson, Dhu, and Schneider (2006), este último generalmente obtenido con Software tales como SHAKE, DEEPSOIL y EERA. Es importante señalar que para simular la propagación de la onda en Comsol Multiphysics, independientemente que el modelo sea horizontalmente infinito este debe ser ampliado incluyendo las regiones de absorción, teoría abordada en la revisión de literatura. La función de transferencia sísmica generada por Comsol y post-procesada en Matlab es una buena aproximación encontrándose en el mismo rango de frecuencia en los distintos modos. A continuación se muestran los 3 modelos más ilustrativos del estudio, los demás modelos pueden ser apreciados en el Anexo 3.

66

Figura 6.3. Modelo 1, Bo. El Bóer: (superior) Funciones de transferencia sísmica generadas por Comsol y Matlab, (centro) perfil 2D y (inferior) Columnas estratigráficas 1D.

Se modelaron los perfiles en 1D con sus respectivas regiones de absorción en Comsol (Figura 6.3 (inferior)) con el fin de obtener las funciones de transferencia sísmica (Figura 6.3 (superior)) las cuales al compararse con la teoría 1D brindan respuesta satisfactoria respecto a los rangos frecuenciales. Por otro lado, las amplificaciones no son completamente similares debido a que Comsol está sustentado en la teoría de elementos finitos siendo un método aproximado dependiente de la discretización del dominio comúnmente conocido como enmallado (mesh); mientras que el algoritmo de respuesta lineal es exacta.

67

Figura 6.4. Modelo 3, Bo. Jorge Dimitrov: (superior) Funciones de transferencia generadas por Comsol y Matlab, (centro) perfil 2D y (inferior) Columnas estratigráficas 1D.

La figura 6.4 (superior) muestra que en contenido de frecuencia ambas respuestas dan prácticamente iguales, antes de la transición trabaja a frecuencias entre 4-6Hz mientras después de la transición se desplaza a frecuencias mayores estando entre los 7-10Hz aproximadamente; esto se sustenta en la teoría que estratos someros tienden a tener altas frecuencias y estratos profundos a frecuencias bajas.

68

Figura 6.5. Modelo 6, Bo. Recreo Sur: (superior) Funciones de transferencia generadas por Comsol y Matlab, (centro) perfil 2D y (inferior) Columnas estratigráficas 1D.

En este modelo se puede apreciar que ambas condiciones (A.T, D.T) se encuentran en rangos de frecuencias similares (4-6Hz). Esto debido a que la diferencia de profundidad entre ambos es mínima (menor a los 5 m) ilustrándose en la figura 6.5 (inferior).

69

6.2.2. Análisis de no-horizontalidad de estratos Para la creación de las funciones de transferencia sísmica hechas en la zona Noroeste de la ciudad de Managua, ya habiendo sido validado el uso de Comsol como una herramienta útil para aproximar los modelos en 1D, se generaron 27 modelos en total (Anexo 1) siendo post-procesados en rutinas de Matlab, de los cuales 18 representan la homogeneidad del suelo en las condiciones (A.T, D.T) y 9 que constituyen la heterogeneidad del suelo (T). Cada línea de transición consta de 3 funciones de transferencias para las condiciones antes mencionadas (A.T, T, D.T) que fueron posteriormente superpuestas. Esto permitió identificar rangos de frecuencias significativos y el posible efecto de la heterogeneidad del suelo al comparar las 3 condiciones. En la Tabla 6.3 se detallan los parámetros necesarios para el modelamientos en Comsol y así obtener las funciones de transferencia sísmica. Tabla 6.3. Parámetros geométricos para modelamiento en Comsol

𝝀 ν (Poisson) Vs (m/s) Vp (m/s) f min (Hertz) (Long. Modelos Onda) Mod. 1 0.30 797.67 1492.30 1 1492.30 Mod. 2 0.30 857.78 1604.75 1 1604.75 Mod. 3 0.30 929.05 1738.09 1 1738.09 Mod. 4 0.30 816.50 1527.53 1 1527.53 Mod. 5 0.30 816.50 1527.53 1 1527.53 Mod. 6 0.30 906.67 1696.21 1 1696.21 Mod. 7 0.30 866.91 1621.83 1 1621.83 Mod. 8 0.30 1013.57 1896.22 1 1896.22 Mod. 9 0.30 923.41 1727.53 1 1727.53

L' (m) H' (m) 4477 4814 5214 4583 4583 5089 4865 5689 5183

7461 8024 8690 7638 7638 8481 8109 9481 8638

En la figura 6.6 se presentan 3 tendencias del resultado de las funciones de transferencia, la primera trata de un comportamiento irregular durante la transición (superior derecho) ubicado en el Bo. El Bóer, el siguiente muestra un comportamiento semi-regular en la transición (centro derecho) situado en el Bo. Jorge Dimitrov y finalmente se visualiza comportamiento aparentemente lineal (inferior derecho) que destina el Bo. El Recreo. Los modelos restantes muestran comportamientos diferenciados entre sí, puesto a su impedancia representativa o por la profundidad somera que muestran.

70

Figura 6.6. Perfiles y Funciones de transferencia sísmica más representativos: (superior) Bo. El Bóer, (centro) Bo. Jorge Dimitrov y (inferior) Bo. El Recreo.

En el modelo 1 correspondiente al Bo. El Bóer, la condición antes de la transición se presentan picos de amplificación en rangos de frecuencia entre 2Hz y 3Hz, lo cual también ocurre en los modelos 4,5 que corresponden a Bo. El Cortijo y M.H.C.P. Para el caso de los modelos 2, 6 ubicados en Bo. William Díaz y Bo. Recreo Sur los rangos de frecuencia significativos oscilan entre 4Hz y 5Hz. Los rangos de 71

frecuencia más altos para esta condición van de 5Hz-7Hz para el modelo 3 y 8 localizados en Bo. Jorge Dimitrov y Bo. La Esperanza, y 8Hz-10Hz para los modelos 7 y 9. En la condición después de la transición los picos de amplificación pueden ser agrupados del siguiente modo, un primer rango entre 8Hz-10Hz para el modelo 1, 3, 5, 8; un segundo rango entre 4Hz-6Hz para el modelo 2, 4, 6 y un tercer rango que oscila entre 2Hz y 3Hz para el modelo 9 en el cual los estratos de suelo son profundos con periodos naturales de vibración altos, y por tanto al ser inversamente proporcional a la frecuencia esta tiende a ser baja. Finalmente un cuarto rango entre 7Hz-9Hz para el modelo 7. Finalmente la función de transferencia sísmica obtenida durante la transición en algunos sitios tiende a ser sumamente irregular, no lográndose identificar picos de amplificación en rangos de frecuencia determinados. Estos son los modelos 1, 2, 5 que se encuentran en El Barrio el Bóer, Barrio William Díaz y M.H.C.P respectivamente, cabe destacar que estos 3 modelos presentan una característica en común, mostrando saltos de profundidad con pendientes mayores al 10% estableciendo las relaciones de estratos más profundos a más someros. Los modelos 4 y 9 que corresponden al Barrio El Cortijo y Planes Altamira, correspondientemente presentan un comportamiento indefinido. En los modelos 6, 7, 8 correspondientes a los sitios: Barrio Recreo Sur, Barrio San José Oriental y Barrio la Esperanza, durante la transición las funciones de transferencia presentan un comportamiento sumamente semejante a las que se obtienen si se asume la condición de horizontalmente infinitos, es importante señalar que a pesar de darse una transición, ésta presenta pendientes menores al 10%, es decir, que la diferencia de niveles de antes y después de la transición es relativamente pequeña (menos de 7m). El modelo 3 localizado en el Barrio Jorge Dimitrov es una especie de eslabón entre los modelos que presentan pendientes mayores al 10% y los que presentan pendientes menores al 10%, el comportamiento de la función de transferencia tiende presentar amplificaciones significativas en un rango de frecuencia especifico entre 7Hz y 9Hz.

72

6.3. ACELEROGRAMAS SINTÉTICOS EN SUPERFICIE Al no disponer de registros de terremotos fuertes en cada punto en estudio, es requerido generarlos por medio de métodos numéricos. La señal sísmica más clara obtenida hasta el momento, data del terremoto de Managua de 1972 registrado por la estación de la refinería ESSO. Debido a que este registro no se obtuvo en la roca elástica es necesario eliminar el efecto del suelo blando haciendo una aproximación para un registro de entrada más o menos libre de amplificaciones. Para este estudio se usa la señal deconvolucionada por Castillo and Zepeda (2013).

Figura 6.7. Acelerograma deconvolucionado del terremoto de Managua de 1972 (Castillo & Zepeda, 2013).

El acelerograma en roca deconvolucionado (Figura 6.7) presenta una aceleración máxima cerca a los 0.23 g, referente a los acelerogramas obtenidos en superficie mediante la convolución, estos presentan máximas aceleraciones en una de cuatro posibles rangos; Un primer rango entre 0.5g y 1g, un segundo rango entre 1g y 1.5g, un tercer rango entre 1.5g y 2g, Finalmente uno aproximado entre 2g y 2.5g. En el Bo. El Bóer (Figura 6.8 superior) las máximas aceleraciones se dan en un tiempo relativamente prolongado (0-10 s), mientras que en el Bo. Jorge Dimitrov (Figura 6.8 centro) el tiempo oscilante del sismo se encuentra entre los 0-8 s; por otro lado en el Bo. Recreo Sur (figura 6.8 inferior) muestra dos comportamientos oscilatorios, primero durante los 0-7 s y luego entre los 10-16 s aproximadamente.

73

Figura 6.8. Acelerogramas sintéticos representativos durante la transición: (superior) Bo. El Bóer, (centro) Bo. Jorge Dimitrov, (inferior) Bo. Recreo Sur.

74

Tabla 6.4. Registro de aceleraciones máximas

MODELO 1

MODELO 3

MODELO 5

MODELO 7

MODELO 9

A.T T D.T A.T T D.T A.T T D.T A.T T D.T A.T T D.T

1g