ESTUDIO NUMÉRICO DE LA PROPAGACIÓN DE...

ESTUDIO NUMÉRICO DE LA PROPAGACIÓN DE ONDASEN SUELOS MEDIANTE FUNCIONES DE GREEN

jorge maestre herediaU

NIVERSITAT

IS

LITTERARIA

EHISPALE

NSIS

SIGILLUM

trabajo fin de máster

dirigida por pedro galvín barrera, antonio romero ordóãsez

doctores por la universidad de sevilla

ÍNDICE GENERAL

Índice de figuras vÍndice de Tablas x1 introducción 1

1.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Estado del conocimiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.4 Organización del texto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 propagación de ondas en el sólido 5

2.1 Ecuaciones generales del movimiento . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Descomposición del vector desplazamiento . . . . . . . . . . . . . . . . . . . 7

2.2.1 Ondas longitudinales y transversales . . . . . . . . . . . . . . . . . . . . . 7

2.2.2 Movimiento en el plano (P-SV) y fuera del plano (SH) . . . . . . . . . . . 8

2.3 Soluciones fundamentales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 solución fundamental del semiespacio 3d homogéneo en el do-minio del tiempo 15

3.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.2 Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.3 Solución en el dominio transformado de Laplace . . . . . . . . . . . . . . . 18

3.3.1 Solución particular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.3.2 Solución homogénea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3.3 Solución en el dominio transformado . . . . . . . . . . . . . . . . . . . . . . 21

3.4 Respuesta sobre la superficie libre . . . . . . . . . . . . . . . . . . . . . . . . 22

3.4.1 Resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.5 Respuesta dentro del semiespacio . . . . . . . . . . . . . . . . . . . . . . . . 42

3.5.1 Resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

4 solución fundamental del semiespacio estratificado en el do-minio del tiempo 57

4.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.2 Solución del semiespacio homogéneo . . . . . . . . . . . . . . . . . . . . . . 59

4.2.1 Problema bidimensional en el plano (SH) . . . . . . . . . . . . . . . . . . . 59

4.2.2 Problema bidimensional anti-plano (P-SV) . . . . . . . . . . . . . . . . . . 63

4.3 Solución de una placa estratificada mediante el TLM . . . . . . . . . . . . . 74

4.3.1 Planteamiento del método en coordenadas cartesianas . . . . . . . . . . . 76

4.3.2 Problema bidimensional anti-plano (SH) . . . . . . . . . . . . . . . . . . . 80

4.3.3 Problema bidimensional en el plano (P-SV) . . . . . . . . . . . . . . . . . . 83

4.4 Solución del semiespacio estratificado . . . . . . . . . . . . . . . . . . . . . . 87

iii

iv índice general

4.4.1 Problema bidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.4.2 Problema tridimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

4.5 Precisión y convergencia de la solución . . . . . . . . . . . . . . . . . . . . . 95

5 ejemplos numéricas 101

5.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

5.2 Validación experimental de la solución fundamental estratificada . . . . . . . 101

5.2.1 Identificación de las propiedades dinámicas del suelo . . . . . . . . . . . 102

5.2.2 Validación numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

5.3 Influencia de la profundidad de la fuente o el receptor . . . . . . . . . . . . 106

5.4 Influencia del amortiguamiento del suelo . . . . . . . . . . . . . . . . . . . . 112

5.5 Suelo compuesto por un estrato blando sobre semiespacio rocoso . . . . . . 115

5.6 Suelo con inclusiones rocosas . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

6 conclusiones 123

6.1 Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Bibliografía 125

ÍNDICE DE F IGURAS

Figura 2.1 Ondas longitudinales y transversales . . . . . . . . . . . . . . . . 9

Figura 3.1 Geometría del semiespacio: superficie libre coincide con el planoz = 0 tal que el eje z positivo define la región del semiespacio.Se evalúan los desplazamientos u en el punto (x, y, z) debido ala fuente b ubicada en la posición (0, 0, z′) . . . . . . . . . . . . . 16

Figura 3.2 Representación de la variable de integración q en el plano com-plejo. Esta representación toma el nombre de Contorno de Cag-niard. Triantafyllidis [48] . . . . . . . . . . . . . . . . . . . . . . . 25

Figura 3.3 Representación del semiespacio para el problema particular enel que receptor se localiza sobre la superficie . . . . . . . . . . . 29

Figura 3.4 Evolución temporal de las componentes de la función Green adi-mensional, GH(2× 103, 0, 0, t; 0, 0, 10 × 103), de un semiespaciohomogéneo. Ángulo de emisión θ = 11.31o. . . . . . . . . . . . . 33

Figura 3.5 Evolución temporal de las componentes de la función Green adi-mensional, GH(10× 103, 0, 0, t; 0, 0, 2 × 103), de un semiespaciohomogéneo. Ángulo de emisión θ = 78.69o. . . . . . . . . . . . . 34

Figura 3.6 Evolución temporal de las componentes de la función Green adi-mensional, GH(10× 103, 0, 0, t; 0, 0, 0.2× 103), de un semiespaciohomogéneo. Ángulo de emisión θ = 88.85o. . . . . . . . . . . . . 35

Figura 3.7 Evolución temporal de las componentes de la función Green adi-mensional, GH(10× 103, 0, 0, t; 0, 0, 1× 10−5), de un semiespaciohomogéneo. Ángulo de emisión θ ≃ 90o. . . . . . . . . . . . . . 36

Figura 3.8 Tipo de singularidad del integrando del término M para ángulode emisión θ = π/2 . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Figura 3.9 Curva de Jordan en el plano complejo p . . . . . . . . . . . . . . 39

Figura 3.10 Comparativa de la función de Green GH33(10× 103, 0, 0, t; 0, 0, 0)

computados mediante el método del aislamiento de la singula-ridad en los extremos de la integral (Método 1), y el métodobasado en la aplicación del Teorema de los residuos (Método 2). . 41

Figura 3.11 Validación de la solución de Johnson con la solución superficialen el plano de Pekeris. a) Evolución temporal del desplazamien-to superficial horizontal. b) Evolución temporal del desplaza-miento superficial vertical . . . . . . . . . . . . . . . . . . . . . . . 41

v

vi Índice de figuras

Figura 3.12 Evolución temporal de la función de Green Gδ11(15, 0, 0, t; 0, 0, 0)

computados por tres formas: Mediante diferencias; medianteintegral de Duhamel para una excitación coseno impulsiva dis-creta; y por el mismo procedimiento continuo . . . . . . . . . . 42

Figura 3.13 Representación del semiespacio del problema general relativo alos términos de las ondas directas P y S . . . . . . . . . . . . . . 43

Figura 3.14 Representación del semiespacio del problema general relativo alos términos de las ondas reflejadas PP y SS . . . . . . . . . . . 45

Figura 3.15 Representación del semiespacio del problema general relativo alos términos de las ondas reflejadas PS y SP . . . . . . . . . . . 47

Figura 3.16 Evaluación de la variable q(p, t) mediante procedimiento numé-rico y analítico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Figura 3.17 Evolución temporal de las componentes de la función Greenadimensional, GH(2× 103, 0, 1× 103, t; 0, 0, 2× 103), de un semi-espacio homogéneo. . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Figura 3.18 Evolución temporal de las componentes de la función Greenadimensional, GH(2× 103, 0, 4× 103, t; 0, 0, 2× 103), de un semi-espacio homogéneo. . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Figura 4.1 Geometría del semiespacio bidimensional: Superficie libre coin-cide con el plano z = 0 tal que el eje z negativo define la regióndel semiespacio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Figura 4.2 Evolución temporal de las componentes de la función Greenadimensional, Gδ(10, 0, 0, t; 0, 0, 0) en el plano (P-SV) de un se-miespacio homogéneo. Linea continua negra corresponde conla función de Green propuesta por Park y Kausel, y los círculosgrises corresponde con la exacta propuesta por Eringen y Suhubi 63

Figura 4.3 Contorno complejo de integración de la función de Green. Parky Kausel [37] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Figura 4.4 Evaluación de los núcleos de las integrales Bxx, Bxz, Bzz para elcaso concreto de z = 0, ν = 0.25. . . . . . . . . . . . . . . . . . . . 71

Figura 4.5 Evaluación de los núcleos de las integrales Bxx, Bxz, Bzz para elcaso concreto de z = 0, ν = 0.25 y obviando el término senoidal(linea continua). En línea discontinua se muestra las asíntotas decada núcleo respectivo. . . . . . . . . . . . . . . . . . . . . . . . . 72

Figura 4.6 Evolución temporal de las componentes de la función Greenadimensional, Gδ(10, 0, 0, t; 0, 0, 0) en el plano (P-SV) de un se-miespacio homogéneo. Linea continua negra corresponde conla función de Green propuesta por Park y Kausel, y los círculosgrises corresponde con la exacta propuesta por Eringen y Suhubi 75

Figura 4.7 Discretización de una lámina. Kausel [22]. . . . . . . . . . . . . . 77

Índice de figuras vii

Figura 4.8 Acoplamiento de n láminas en el sentido de los elementos fini-tos. Kausel [24]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

Figura 4.9 Evolución temporal de las componentes de la función Greenadimensional, Gyy

(5, 0, 0, t; 0, 0, 0) en el plano (SH) de una placahomogénea. En línea continua se dibuja la solución del TLM yen círculos la solución exacta . . . . . . . . . . . . . . . . . . . . 83

Figura 4.10 Evolución temporal de las componentes de la función Greenadimensional, Gzz(5, 0, 0, t; 0, 0, 0) en el plano (P-SV) de una pla-ca homogénea. En línea continua se dibuja la solución del TLMy en círculos la solución exacta . . . . . . . . . . . . . . . . . . . 87

Figura 4.11 Superposición de la respuesta de la placa estratifica con el semi-espacio homogéneo. Kausel [24]. . . . . . . . . . . . . . . . . . . 88

Figura 4.12 Estabilidad del método híbrido para tres pasos temporales ∆t =1/(2 fmax), ∆t = 1/(4 fmax) y ∆t = 1/(8 fmax). Función de Greendel semiespacio ante una carga vertical superficial tipo impulsi-va senoidal. Gzz(1000, 0, 0, t, 0, 0, 0). En línea continua se dibujala respuesta del método híbrido y en círculos grises la respuestaexacta. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

Figura 4.13 Precisión del método híbrido para tres discretizaciones Nλ =

2,Nλ = 4 y Nλ = 8. Función de Green del semiespacio ante unacarga vertical superficial tipo impulsiva senoidal. Gzz(1000, 0, 0, t, 0, 0, 0).En línea continua se dibuja la respuesta del método híbrido yen círculos grises la respuesta exacta. . . . . . . . . . . . . . . . . 99

Figura 5.1 Espectro de onda en el dominio número de onda-frecuencia pa-ra el ensayo experimental y para los modelos teóricos. . . . . . 105

Figura 5.2 Comparación de las curvas de dispersión experimental (líneanegra), con las curvas teóricas del modelo 1 (linea gris claro) ydel modelo 2 (línea gris oscura) . . . . . . . . . . . . . . . . . . . 106

Figura 5.3 Comparación de las curvas de dispersión experimental (líneanegra), con las curvas teóricas del modelo 2 para las distanciasde medida máxima de xmax = 84m (linea gris oscuro), xmax =66m (linea gris claro) y xmax = 48m (linea gris discontinua) . . . 107

viii Índice de figuras

Figura 5.4 Comparativa entre la solución del semiespacio y la solución delespacio completo. La fuente y el receptor se encuentran a lamisma cota y a 5 m de distancia. Se representa nueve curvasrelativas a las profundidades: d = 0 m, d = 1 m, d = 2 m, d =4 m, d = 8 m, d = 16 m, d = 32 m, d = 64 m, d = 128 m.En linea negra continua se representa la evolución temporal delos desplazamientos calculados con la expresión de Johnson, enlinea discontinua gris oscuro la correspondiente al espacio com-pleto, y en gris claro discontinuo con marcas se representa eltiempo de llegada de las ondas. . . . . . . . . . . . . . . . . . . . 109

Figura 5.5 Comparativa entre la solución del semiespacio y la solución delespacio completo. La fuente y el receptor se encuentran a lamisma cota y a 50 m de distancia. Se representa nueve curvasrelativas a las profundidades: d = 0 m, d = 1 m, d = 2 m, d =

4 m, d = 8 m, d = 16 m, d = 32 m, d = 64 m, d = 128 m.En linea negra continua se representa la evolución temporal delos desplazamientos calculados con la expresión de Johnson, enlinea discontinua gris oscuro la correspondiente al espacio com-pleto, y en gris claro discontinuo con marcas se representa eltiempo de llegada de las ondas. . . . . . . . . . . . . . . . . . . . 110

Figura 5.6 Comparativa entre la solución del semiespacio y la solución delespacio completo. La fuente y el receptor se encuentran a lamisma cota y a 100 m de distancia. Se representa nueve curvasrelativas a las profundidades: d = 0 m, d = 1 m, d = 2 m, d =4 m, d = 8 m, d = 16 m, d = 32 m, d = 64 m, d = 128 m.En linea negra continua se representa la evolución temporal delos desplazamientos calculados con la expresión de Johnson, enlinea discontinua gris oscuro la correspondiente al espacio com-pleto, y en gris claro discontinuo con marcas se representa eltiempo de llegada de las ondas. . . . . . . . . . . . . . . . . . . . . 111

Figura 5.7 Evolución de los desplazamientos verticales debido a una exci-tación impulsiva vertical para diferentes factores de amortigua-miento. La fuente y el receptor se encuentra sobre la superficiey a 6 m de distancia. Se representa once curvas relativas a losfactores de amortiguamiento: ξ = 0, ξ = 0.01, ξ = 0.02, ξ =

0.03, ξ = 0.04, ξ = 0.05, ξ = 0.06, ξ = 0.07, ξ = 0.08, ξ =0.09, ξ = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Figura 5.8 Evolución de los desplazamientos verticales debido a una excita-ción tipo escalón vertical para diferentes factores de amortigua-miento. La fuente y el receptor se encuentra sobre la superficiey a 6 m de distancia. Se representa once curvas relativas a losfactores de amortiguamiento: ξ = 0, ξ = 0.01, ξ = 0.02, ξ =

0.03, ξ = 0.04, ξ = 0.05, ξ = 0.06, ξ = 0.07, ξ = 0.08, ξ =0.09, ξ = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figura 5.9 Evolución de los desplazamientos verticales con respecto al tiem-po (eje x) y en función de la distancia a la fuente (eje y) debido auna carga impulsiva. La fuente y el receptor se encuentra sobrela superficie. Se ha considerado un amortiguamiento materialnulo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

Figura 5.10 Evolución temporal de los desplazamientos verticales de un se-miespacio estratifica compuesto por un estrato blando de 5 mde espesor sobre un suelo rocoso debido a una excitación im-pulsiva tipo senoidal vertical. La fuente y el receptor se en-cuentra sobre la superficie y separados una distancia variable(x = 2, 4, 6, 8, 10 m). . . . . . . . . . . . . . . . . . . . . . . . . . 117

Figura 5.11 Comparativa de la evolución temporal de los desplazamientosde un semiespacio estratificado (linea continua) respecto a losdesplazamiento de un estrato libre-libre (linea discontinua gris).El semiespacio está compuesto por un estrato blando de 5 mde espesor sobre un suelo rocoso, mientras que la placa libretiene el mismo espesor y propiedades que el estrato blando. Laexcitación es superficial vertical impulsiva de tipo seniodal y seregistra los desplazamientos verticales superficiales a 4 m dedistancia. En linea de punto y raya gris claro se representa lainversa de los desplazamientos debido a las reflexiones en laplaca. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

Figura 5.12 Evolución temporal de los desplazamientos verticales de un se-miespacio estratifica debido a una excitación impulsiva tipo se-noidal vertical. El semiespacio está compuesto suelo blando enel que se le intercala un delgado estrato rocoso de 0.5 m deespesor y a 5 m de profundidad. La fuente y el receptor se en-cuentra sobre la superficie y separados una distancia variable(x = 2, 4, 6, 8, 10 m). . . . . . . . . . . . . . . . . . . . . . . . . 120

ix

Figura 5.13 Comparativa de la evolución temporal de los desplazamientosde un semiespacio estratificado compuesto por un suelo blandoen el que se le intercala un estrato rocoso de 0.5 m de espesory a 5 m de profundidad (linea continua), respecto de un sueloblando de 0.5 m de espesor que apoya sobre semiespacio ro-coso infinito. Las propiedades de los estratos son idénticas. Laexcitación es superficial vertical impulsiva de tipo seniodal y seregistra los desplazamientos verticales superficiales a 4 m dedistancia. En linea de punto y raya gris claro se representa lainversa de los desplazamientos debido a las reflexiones en laplaca. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

ÍNDICE DE TABLAS

Tabla 4.1 Coste adicional computacional. . . . . . . . . . . . . . . . . . . . 72

Tabla 4.2 Error ramas integrales . . . . . . . . . . . . . . . . . . . . . . . . 73

Tabla 4.3 Funciones asintóticas . . . . . . . . . . . . . . . . . . . . . . . . . 73

Tabla 4.4 Matrices elementales del TLM. Problema plano SH. . . . . . . . . 81

Tabla 4.5 Matrices elementales del TLM. Problema plano P-SV. . . . . . . 84

Tabla 5.1 Perfil estratigráfico . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Tabla 5.2 Perfil estratigráfico . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Tabla 5.3 Tipos de ondas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

Tabla 5.4 Tipos de ondas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

ACRÓNIMOS

La forma plural se construirá añadiendo una s al final del acrónimo.

x

1INTRODUCCIÓN

1.1 introducción

El presente trabajo pretende implementar una herramienta fiable y eficiente con laque poder predecir las vibraciones producidas por cargas dinámicas. Este problema sepuede dividir en tres bloques. En el primero se debe analizar y modelizar los mecanis-mos de emisión de vibraciones como por ejemplo las vibraciones causadas por el pasodel tren, o las ocasionadas por explosivos durante la ejecución de túneles. En segundolugar se requiere estudiar la propagación de ondas en el suelo. La propagación de lasondas depende de las propiedades dinámica del suelo, de su estratificación y de lascondiciones de contorno. Finalmente estas vibraciones, que suponen un campo de on-das incidentes, pueden afectar a las estructuras, al funcionamiento de sus instalacioneso a sus ocupantes, lo que requiere estudiar la interacción suelo-estructura.

Este documento se centra en la parte de propagación de ondas en el suelo, dondelos modelos numéricos basados en la formulación del Método de los Elementos deContorno (MEC) en el dominio del tiempo permite estudiar con rigor y eficiencia estetipo de problema [15, 29, 30]. El MEC se basa en una ecuación integral de superficie,y relaciona un estado conocido virtual, que corresponde con la solución fundamental,con uno desconocido real. La diferencia y ventaja respecto otros método numéricos dedominio, como el Método de los Elementos Finitos (MEC), es que solo requiere discre-tizar el contorno, en concreto la parte del mismo en el que las condiciones de contornodel problema en cuestión son diferentes a las de la solución fundamental. Ademáspara medios infinitos o semi-infinitos satisface implícitamente las condiciones de ra-diación. Sin embargo, una de las limitaciones del método es que no existen solucionesfundamentales para todos los problemas, o a veces la solución no se encuentra enforma explícita. De hecho es común utilizar la solución fundamental del semiespaciocompleto para problemas elastodinámicos, debido a que existe en forma explícita yes sencillo introducirlo en la formulación del MEC. Por contrapartida, para resolverproblemas de propagación de ondas en el suelo se requiere mallar la superficie e im-poner la condición de tensión nula. En determinados problemas, ello puede suponerun fuerte impedimento por el alto coste computacional que conlleva y a veces haceque sea inviable directamente su aplicación. En este campo es interesante la soluciónfundamental del semiespacio, en donde la superficie no requiere ser mallada. Aquí esdonde nace el interés del presente trabajo que consiste en estudiar la solución funda-mental del semiespacio homogéneo y estratificado, implementarla y validarla, con elobjetivo de que pueda ser introducida posteriormente en la formulación del MEC.

1

2 introducción

1.2 estado del conocimiento

Desde que Lamb [27] en 1904 planteara el problema de una carga concentrada eimpulsiva sobre un semiespacio homogéneo, isótropo, elástico y lineal, numeroso in-vestigadores han tratado de encontrar la solución. Comenzando por Pekeris [38] en1955 y Chao [9] en 1960 obtuvieron de forma explícita la respuesta superficial anteuna carga puntual de tipo escalón vertical y horizontal sobre la superficie de un semi-espacio homogéneo con coeficiente de Poisson de 0.25. Más tarde en 1974 Mooney [31]generalizó la resultados de Pekeris para un coeficiente de Poisson arbitrario. En 1975

Eringen and Suhubi [16] presentaron la respuesta en el interior del semiespacio parauna carga impulsiva superficial vertical y horizontal para el caso bidimensional, y parael caso tridimensional la respuesta en la superficie hasta un coeficiente de Poisson de0.26. En la misma época, Johnson [19] estudió el caso general para cualquier coeficien-te de Poisson ,ubicación y dirección tanto de la fuente como del receptor, y obtuvola respuesta en forma integral. Actualmente, esta última es la única solución generaltridimensional analítica en el dominio temporal que se conoce. Aunque es cierto queexisten soluciones generales en el dominio transformado número de onda-frecuencia[26].Por otro lado, para el problema del semiespacio horizontalmente estratificado, no

existen soluciones directamente en el dominio espacio-tiempo, sino las soluciones seplantean en el dominio transformado número de onda-frecuencia. En 1981 Bouchon [7]propuso un método original en el que la respuesta tridimensional ante una carga con-centrada impulsiva se puede aproximar por la superposición de ondas producidas poruna cierta distribución de fuentes virtuales alrededor de la fuente original. Físicamentelas fuentes virtuales representan las reflexiones de las ondas sobre las superficies de losestratos. Esta solución solo es válida para una franja de tiempo determinado. En 1983

R. Aspel y E. Luco [3, 4] presentaron la solución al problema axisimétrico medianteel planteamiento de un sistema de ecuaciones que compatibiliza los desplazamientosy tracciones en la interfase de los estratos. La solución de cada estrato se obtiene dela ecuación de Navier transformada al dominio número de onda-frecuencia mediantela trasformada de Hankel y de Fourier. Paralelamente, Kausel y Roësset [26] propu-sieron un método potente y sistemático basado en la matriz de rigidez para resolverel problema tridimensional ante una fuente genérica ubicada en cualquier parte delsemiespacio. El método parte de calcular las soluciones analíticas a la ecuación de Na-vier en el dominio transformado (transformada de Hankel y Fourier) para el problemadel semiespacio y del estrato homogéneo tanto en desplazamientos como en traccio-nes en las superficies libres. Estas soluciones convenientemente dispuestas permitenformular la matriz de rigidez elemental, y expresan la relación entre los desplazamien-tos y las tracciones de las interfases. La respuesta del semiespacio horizontalmenteestratificado se realiza tras ensamblar las matrices elementales y formar una matriz derigidez global. Este método posee cierto parecido al Método de los Elementos Finitos

1.3 objetivos 3

pero la solución lograda en principio es exacta. Mediante esta técnica puede resolver-se una amplia variedad de problemas dinámicos. La respuesta que se obtiene es enel dominio transformado número de onda-frecuencia, lo que requiere una integraciónsobre el número de onda y sobre la frecuencia que numéricamente conlleva ciertosproblemas y necesita de la práctica de determinadas técnicas que aproximan la solu-ción. Como alternativa, posteriormente en 1994 Kausel [22] propuso un nuevo método,llamado Thin-Layer Method (TLM), para obtener las funciones de Green directamenteen el dominio temporal, que es especialmente adecuado para problemas dinámicosresueltos mediante el MEC en el dominio del tiempo. Esta técnica es similar al métododirecto de la matriz de rigidez pero a diferencia de éste la matriz de rigidez provienede una discretización vertical del semiespacio estratificado, donde los desplazamien-tos se aproximan mediante funciones de forma. La matriz a la que se llega es muchomás sencilla matemáticamente, y mediante la técnica de la superposición modal enfrecuencia se obtiene la función de Green tras resolver un problema estándar de auto-valores y autovectores. De esta forma la respuesta se presenta en el dominio númerode onda y tiempo, tal que solo requiere de una integración sobre el número de onda.En principio esta técnica solo es válida para una composición de estratos bien libreso bien apoyados sobre un superficie infinitamente rígida. Sin embargo puede usarseen combinación con la solución del semiespacio homogéneo en el dominio temporaldando lugar a un método híbrido propuesto por Park and Kausel [24, 25].

1.3 objetivos

El propósito de este trabajo es implementar la solución fundamental del semiespacioque luego será introducido en un código donde se usa el MEC el FEM para resolverproblemas dinámicos de interacción suelo-estructura. En este trabajo se presenta laformulación de Johnson para el semiespacio homogéneo en donde las expresiones semuestran en forma integral. Conscientes que el MEC requiere un uso repetitivo de lasolución fundamental, en este trabajo se ha puesto especial atención en la eficienciadel código para minimizar el tiempo de computo y evitar problemas numéricos. Porotro lado, atendiendo a la habitual heterogeneidad del suelo, se ha implementadola solución fundamental propuesta por Kausel para un semiespacio estratificado. Sebasa en un método híbrido que combina la respuesta de una placa estratifica con la delsemiespacio homogéneo. La formulación requiere de una discretización espacial en elsentido de la estratificación y una discretización temporal. Ambas deben ser elegidascon cuidado para que el método sea estable.

1.4 organización del texto

La organización del texto es la siguiente:

4 introducción

Capítulo 1: Introducción. Se realiza una introducción del trabajo, se presenta losobjetivos y se organiza el texto

Capítulo 2: Propagación de ondas en el sólido. Se expone la base matemática sobrepropagación de ondas en un medio isótropo, homogéneo, elástico y lineal.

Capítulo 3: Solución fundamental del semiespacio 3D homogéneo en el dominio del tiem-

po. Se presenta la formulación de Johnson [19] del semiespacio homogéneo. Lasexpresiones se encuentran en forma integral con numerosos singularidades yproblemas numéricos, sobre todo cuando la fuente y el receptor tienden a la su-perficie. En este capítulo se exponen métodos para evitar las singularidades ypaliar el efecto de las ondas superficiales. La implementación de la formulaciónse valida con la solución de Pekeris [38].

Capítulo 4: Solución fundamental del semiespacio estratificado en el dominio del tiempo.Se muestra la formulación propuesta por Park y kausel [24, 25] para un semies-pacio estratificado horizontalmente. Se basa en un método híbrido que combinala solución de una placa horizontalmente estratificada con la solución del semi-espacio en el dominio número de onda-tiempo. Ambos soluciones se acoplanresolviendo un sistema de ecuaciones recursivo para cada pasa de tiempo.

Capítulo 5: Ejemplos numéricas. Una vez presentada la solución fundamental delsemiespacio estratificado se valida con los resultados experimentales, y se resuel-ven numéricamente algunos ejemplos de interés.

Capítulo 6: Conclusiones y trabajo futuro. Por último se resumen las conclusionesdel presente trabajo y se muestra la futura aplicación de las soluciones imple-mentadas.

2PROPAGACIÓN DE ONDAS EN EL SÓL IDO

El comportamiento de un suelo ante cargas dinámicas suele ser, en términos gene-rales, complejo. El suelo está compuesto por una estructura sólida (formada normal-mente por diferentes tipos de materiales) porosa ocupada total o parcialmente poragua. Ello promueve que el comportamiento general del suelo sea anisótropo y nolineal, presentando para el caso de suelos saturados o parcialmente saturados, un aco-plamiento entre la estructura sólida y el agua [5, 6]. No obstante para suelos secos nocohesivos el comportamiento no lineal puede ser despreciado, y puede asumirse uncomportamiento elástico lineal.

En este capítulo, se presenta la base matemática sobre la propagación de ondasen un suelo seco no cohesivo. Se asumen las hipótesis de pequeñas deformaciones ycomportamiento elástico lineal. Haciendo un tratamiento adecuado de las ecuacionesse demuestra que el movimiento de las ondas puede ser descompuesto en un movi-miento longitudinal y otro transversal al frente de onda. Éste último, para el caso quese estudia de un semiespacio homogéneo u horizontalmente estratificado, se puededescomponer a su vez en un movimiento transversal paralelo y normal a la superfi-cie libre. Este hecho tiene relativa importancia puesto que desacopla las ecuaciones ysimplifica su resolución.

2.1 ecuaciones generales del movimiento

Tal como se ha comentado inicialmente, se restringe el campo de estudio a suelos enlos que se supone un comportamiento homogéneo, elástico y lineal. Esta suposicióncabe a priori pensar que es muy restrictiva, si bien puede aplicarse a una gran variedadde suelos. Además la heterogeneidad puede ser tenida en cuenta mediante la uniónde diferentes estratos homogéneos.

Las ecuaciones del movimiento se obtienen por simple aplicación de la teoría de laelasticidad bajo las hipótesis de pequeñas deformaciones y desplazamiento [33, 13]. Seplantean en coordenadas cartesianas y en notación indicial. No obstante, pueden serfácilmente extrapoladas a coordenadas cilíndricas o esféricas.

El tensor de deformaciones de Green para pequeñas deformaciones se expresa enfunción de los desplazamientos ui de la siguiente manera:

ǫij =12(ui,j + uj,i) (2.1)

Esta ecuación es también conocida por ecuación de compatibilidad donde la comadenota derivación espacial.

5

6 propagación de ondas en el sólido

Por otro lado, para un material lineal elástico isótropo, el tensor de tensiones deCauchy puede expresarse en función de las deformaciones de la siguiente forma:

σij = 2µǫij + λǫkkδij (2.2)

que se conoce como ecuación de Lamé o ecuación de comportamiento, donde µ y λ

son los coeficiente Lamé.Para cerrar el problema falta por plantear la ecuación de equilibrio:

σij,j + ρbi = ρui (2.3)

En lo que se refiere a problemas de propagación de ondas, es usual expresar lasecuaciones en función del vector desplazamiento, que es la variable de interés. Paraello basta con sustituir la Ecuación de compatibilidad (2.1) en la Ecuación de com-portamiento (2.2) y a su vez ésta en la Ecuación de equilibrio (2.3) obteniéndose lasEcuaciones de Navier:

µui,jj + (λ + µ)uj,ij + bi = ρui (2.4)

que representa un conjunto de 3 ecuaciones donde el subíndice i es libre. Esta ecuaciónpuede expresarse en notación de teoría de campos de la forma:

µ∇2u+ (λ + µ)∇∇ · u+ b = ρu (2.5)

Sin embargo suele utilizarse una variante de esta ecuación (2.5) que suele simplificarla resolución de los problemas elastodinámicos:

(λ + 2µ)∇∇ · u− µ∇×∇× u+ b = ρu (2.6)

La solución de las ecuaciones de Navier produce un campo de desplazamientos queincluye una serie de constantes, y que hay que determinar aplicando las condicionesiniciales y condiciones de contorno que particularizan el problema en cuestión.Las condiciones iniciales definen el movimiento y la velocidad del sólido para el

instante inicial:

ui(x, 0) = ui(x)

ui(x, 0) = ˙ui(x) (2.7)

Mientras que las condiciones de contorno restringen el campo de desplazamientos ytensiones en el contorno Γ. Si se divide éste en dos en función de que se conozcan lastensiones (Γt) o los desplazamientos (Γu), las condiciones quedan como:

ui(x, t) = ui(x, t) ∀x ∈ Γu

ti(x, t) = ti(x, t) ∀x ∈ Γt (2.8)

donde el vector ti representa las tracciones en el contorno de normal exterior ni, y quepuede calcularse mediante el principio de Cauchy como ti = σijnj.

2.2 descomposición del vector desplazamiento 7

Usando las ecuaciones de Navier junto con las condiciones iniciales y de contorno,el problema queda cerrado. Sin embargo, estas ecuaciones se presentan en forma dife-rencial, lo que origina en sí un problema esencial a la hora de obtener analíticamenteel campo de desplazamientos por integración de las ecuaciones de movimiento enproblemas elasto-dinámicos. En este sentido, Poisson [40, 41] fue el precursor al intro-ducir en 1829 el primer teorema general para la integración de las ecuaciones. Desdeentonces, numerosos investigadores han presentado diferentes técnicas matemáticasque transforman el sistema de ecuaciones de Navier y simplifican su resolución. Caberesaltar la técnica de los potenciales introducidas por Lamé [1] en 1852 que es una delas más simples y útiles. Se basa en la idea de la descomposición de Helmholtz delvector desplazamiento (y vector de fuerza volumétrica) en función de potenciales, deforma que las ecuaciones de Navier quedan desacopladas en un sistema de dos ecua-ciones en función los mismos. Contemporánea a ésta y análoga también debe citarseel trabajo desarrollado por Stokes [46] en 1849 que plantea las ecuaciones de Navier entérminos del movimiento dilatacional (θ = ∇ · u) y rotacional (! = ∇× u), llegándosea un sistema similar. Posteriormente surgieron otras técnicas basadas en la transfor-mación de Fourier del dominio natural espacio-tiempo a número de onda-frecuencia[7, 3, 4], o el método de la transformada de Laplace [11, 19]. Estas técnicas se explicaráncon mayor detalles a lo largo del presente documento.

Si bien, con las ecuaciones mostradas el problema está cerrado, para el caso dedominios infinitos (semiespacio o espacio completo) existe una restricción adicionalimplícita al sistema de ecuaciones [14]. Esta restricción es importante y condicionael comportamiento de los desplazamientos y tracciones en el infinito. Supóngase uncontorno ficticio esférico de radio R suficientemente grande, si se desacopla los des-plazamientos y tracciones en una parte irrotacional y equivolumial (superíndice p y srespectivamente) se tienen la siguiente relaciones:

tpi + ρCpu

pi = O(

1R)

tsi + ρCsusi = O(

1R) (2.9)

que se conoce como condición de radiación donde Cp y Cs son respectivamente lavelocidad de propagación de las ondas irrotacionales y equivolumiales. Esta condiciónobliga a que el flujo de energía fluya hacia el infinito tendiendo a cero asintóticamente.

2.2 descomposición del vector desplazamiento

2.2.1 Ondas longitudinales y transversales

En la mayoría de los problemas elastodinámicos el nivel de dificultad a la hora deobtener una solución analítica suele ser alto e incluso a veces inabordable. Asumiendonulo las fuerzas de volumen (caso habitual), la forma clásica de usar el sistema de


ecuaciones de Navier consiste en desacoplarlo en ecuaciones de Helmholtz medianteel método de los potenciales, lo que reduce considerablemente la complejidad. Éstapráctica es bien utilizada por Kausel [23] y Domínguez [14] que son las referenciasprincipales para este capítulo. Consiste en descomponer el vector desplazamiento endos términos: uno que deriva del gradiente de una función escalar, y otro que provienedel rotacional de una función vectorial, que en notación vectorial se expresa:

u = ∇Φ +∇× Ψ (2.10)

Esta descomposición matemática tiene un significado físico: el primer término hacereferencia al movimiento dilatacional de las partículas, mientras que el segundo almovimiento equivolumial.Aplicando la descomposición (2.10) al sistema de Ecuaciones de Navier (2.6), resulta

el siguiente sistema desacoplado:

(λ + 2µ)∇2Φ = ρΦ (2.11)

µ∇2Ψ = ρΨ (2.12)

que constituyen una ecuación escalar y vectorial respectivamente. De acuerdo conAchenbach [2] a ésta última ecuación hay que añadirle la condición ∇ · Ψ = 0 querelaciona las tres funciones potenciales del vector (Ψx, Ψy y Ψz).Por tanto, para un movimiento armónico se llega a un sistema de dos ecuaciones de

Helmholtz:

∇2Φ + k2PΦ = 0, kP = ω/CP Cp =√

(λ + 2µ)/ρ (2.13)

∇2Ψ + k2SΨ = 0, kS = ω/CS Cp =

√

µ/ρ (2.14)

Desde un punto de vista físico, la Ecuación (2.13) o (2.11) describe el movimientode las ondas de dilatación, y corresponden con un desplazamiento longitudinal delas partículas en el sentido de la dirección del frente de onda. Se conocen tambiéncomo ondas primarias (puesto que son las más rápidas propagándose a velocidad CP)y están asociadas a un cambio de volumen. Por otro lado, la Ecuación (2.12) o (2.14)relativa a las ondas equivolumiales, se caracterizan por un movimiento transversal decizalladura y están asociadas a un cambio de forma. Estas ondas son más lentas, deahí la terminología de ondas secundarias, y se propagan a velocidad CS. Véase figura(2.1)

2.2.2 Movimiento en el plano (P-SV) y fuera del plano (SH)

El movimiento transversal de las partículas puede descomponerse a su vez en dosondas relativas a las dos direcciones del plano transversal [14]. En este punto, se debe


ONDAS LONGITUDINALES

ONDAS TRANSVERSALES

ONDAS LONGITUDINALES

ONDAS TRANSVERSALES

Figura 2.1: Ondas longitudinales y transversales

remarcar que el presente trabajo está orientado a obtener la respuesta de un semi-espacio homogéneo u horizontalmente estratificado. Por lo que asumiendo el planohorizontal de interfase como referencia (x-z), las ondas transversales pueden descom-ponerse en una componente paralela al plano y otra normal, que reciben el nombrede ondas horizontal (SH) y verticalmente porlarizadas (SV) respectivamente. Esta di-visión tiene su sentido puesto que guardan relación con las funciones potenciales an-teriormente introducidas.

Asumiendo el eje z perpendicular al plano horizontal de interfase, de acuerdo a Pi-lant [39] el vector desplazamiento admite una descomposición alternativa de la forma:

u = ∇Φ − l∇×∇× (ezΨ′) +∇× (ezχ) (2.15)

donde Φ Ψ′ χ son funciones escalares o potenciales, y l es un factor que hace coheren-te las dimensiones entre los diferentes potenciales. Introduciendo en la Ecuación deNavier (2.6), se llega nuevamente a un sistema desacoplado:

(λ + 2µ)∇2Φ = ρΦ (2.16)

µ∇2Ψ′ = ρΨ′ (2.17)

µ∇2Ø = ρχ (2.18)

En los que las Ecuaciones (2.16) (2.17) (2.18) hacen referencias las ondas P, SV y SHrespectivamente.

Se restringe el análisis al caso de deformación plana y se estudia el problema depropagación bidimensional en el plano (x-z). Es obvio que en este caso las deriva-


das respecto a la dirección y son nulas. Mediante algunas operaciones algebraicas yrenombrando el potencial como Ψ = l∂Ψ′/∂x la Ecuación (2.15) se reduce a:

u = ∇Φ +∇× (eyΨ) +∇× (ezχ) (2.19)

o lo que es lo mismo:

u =

[∂Φ

∂x− ∂Ψ

∂z

]

~ex

− ∂χ

∂x~ey

+

[∂Φ

∂z− ∂Ψ

∂x

]

~ez (2.20)

que conlleva a que el movimiento en el plano (ux y uz) se desacople completamentedel movimiento fuera del plano (uy). Esto es, el movimiento en el plano queda descritoen función de los potenciales Φ y Ψ que corresponden con las ondas de presión y ver-ticalmente polarizadas (P-SV) respectivamente, y el antiplano en función del potencialχ asociado a las ondas horizontalmente polarizadas. Esto último sólo ocurre en el casobidimensinal.Si introducimos nuevamente la descomposición (2.19) en la Ecuación de Navier (2.6)

se llega finalmente al mismo sistema de ecuaciones pero en función del potencial Ψ:

(λ + 2µ)∇2Φ = ρΦ (2.21)

µ∇2Ψ = ρΨ (2.22)

µ∇2Ø = ρχ (2.23)

Curiosamente, la solución para el caso tridimensional en coordenadas cilíndricasde un semiespacio homogéneo u horizontalmente estratificado es totalmente análoga,donde ahora el vector desplazamiento se formula según las componentes radial, an-gular y vertical (ur, uθ, uz). El desarrollo de la solución se realiza siguiendo la mismalínea que para el problema bidimensional, e igualmente se llega a un sistema de tresecuaciones escalares desacoplados en función de los potenciales Φ,Ψ,χ asociados alas ondas P, SH y SV respectivamente. En este caso las componentes del vector despla-zamiento no quedan desacopladas, dependiendo cada una de los tres potenciales. Noobstante los desplazamientos pueden obtenerse mediante combinación adecuadas delas soluciones en el plano y fuera del plano. Un desarrollo más explicito del problematridimensional puede encontrarse en [23] y [43]

2.2.2.1 Ondas de presión

Tal como se ha deducido el movimiento de propagación se formula mediante un sis-tema de ecuaciones escalares desacoplados. Cada ecuación puede resolverse indepen-dientemente, de modo que el vector desplazamiento se expresa como una combinación


de las soluciones de los potenciales. Sólo para el caso bidimensional esta combinaciónse encuentra desacoplada por un vector desplazamiento en el plano y fuera del plano.

En lo que se refiere a la Ecuación (2.21) describe el movimiento longitudinal de lasparticulas en el plano. Esta ecuación puede expresarse de la forma de una ecuaciónde Helmholtz, y su resolución puede abordarse de diferentes maneras. A continua-ción se presenta el método de la transformada de Fourier, cuyo método es utilizadoextensamente por Kausel [23], de forma que la ecuación se resuelve en el dominiosimplificado número del número de onda-frecuencia, y luego se devuelve al dominionatural espacio-tiempo mediante una serie de antitransformadas.

Realizando una transformada del tiempo a frecuencia la ecuación (2.21) se reduce a:

(λ + 2µ)∇2Φ + ω2ρΦ = 0 (2.24)

o abreviadamente:

∇2Φ + kPΦ = 0 kP = ω/CP Cp =√

(λ + 2µ)/ρ (2.25)

donde se ha eliminado el diferencial respecto del tiempo. El sombrero sobre la variableindica que se trabaja en el dominio de la frecuencia.

Asumiendo ahora que la geometría es invariable en dirección x, tras realizar unanueva transformación del espacio al número de onda, la ecuación se simplifica a:

k2xΦ +d2Φ

dz2+ kpΦ = 0 (2.26)

donde la tilde sobre la variable significa que se trabaja en el dominio número de onda-frecuencia.

Por tanto se llega a una ecuación diferencial ordinaria de segundo orden en direc-ción z cuya solución genérica es de la forma:

Φ(kx, z,ω) = Ipe−ikzpz + Rpe

+ikzpz (2.27)

Introduciendo la solución (2.27) en la Ecuación (2.26) se deduce la relación entre el nú-mero de onda vertical (kzp) y el número de onda horizontal (kx) para cada frecuencia:

k2zp + k2x = k2p =ω2

C2P

(2.28)

Las amplitudes IP y Rp son constantes de integración y físicamente hacen referenciala parte incidente (que se propagan en dirección positiva del eje z) y reflectadas (quese propagan en dirección negativa al eje z) de las ondas P. Para Kpz siendo real corres-ponde con la propagación de la onda en dirección ~ez, mientras que siendo imaginariasignifica la evanescencia de la onda que decae exponencialmente.

La solución del potencial en el dominio natural se obtiene sin más que aplicar laantitransformada de Fourier dos veces:

Φ(x, z, t) =1

(2π)2

∫ ∞

−∞

∫ ∞

−∞ei(ωt−kxx)

(

Ipe−ikzpz + Rpe

+ikzpz)

dkxdω (2.29)


2.2.2.2 Ondas verticalmente polarizadas

La Ecuación (2.22), que hace referencia al movimiento vertical en el plano, se puederesolver de forma totalmente similar al apartado anterior. Realizando un doble trans-formación al dominio del número de onda y de la frecuencia se tiene:

k2xΨ +d2Ψ

dz2+ ksΨ = 0 (2.30)

Cuya solución es:Φ(kx, z,ω) = ISVe

−ikzsz + RSVe+ikzsz (2.31)

con

k2zs + k2x = k2s =ω2

C2S

(2.32)

Las amplitudes ISV y RSV son constantes de integración y tienen una interpretaciónanáloga, esto es, se refieren a la parte incidente y reflectadas de las ondas SV.La solución en el dominio natural se obtiene al aplicar la anti-transformada de Fou-

rier a la Ecuación (2.31):

Ψ(x, z, t) =1

(2π)2

∫ ∞

−∞

∫ ∞

−∞ei(ωt−kxx)

(

ISVe−ikzsz + RSVe

+ikzsz)

dkxdω (2.33)

2.2.2.3 Ondas horizontalmente polarizadas

El movimiento fuera del plano de las ondas horizontalmente polarizadas viene go-bernado por al ecuación (2.23). Tal como se ha practicado anteriormente se aborda suresolución mediante una doble transformación, llegándose en este caso a la ecuación:

k2xχ +d2χ

dz2+ ksχ = 0 (2.34)

Cuya solución se expresa de la forma:

χ(kx, z,ω) = ISHe−ikzsz + RSHe

+ikzsz (2.35)

donde se tiene la siguiente relación del número de onda vertical y horizontal:

k2zs + k2x = k2s =ω2

C2S

(2.36)

Las amplitudes ISH y RSH son constantes de integración y se refieren a la parte inci-dente y reflectadas de las ondas SH.La solución en el dominio natural se obtiene al aplicar una doble anti-transformada

de Fourier a la Ecuación (2.35):

χ(x, z, t) =1

(2π)2

∫ ∞

−∞

∫ ∞

−∞ei(ωt−kxx)

(

ISHe−ikzsz + RSHe

+ikzsz)

dkxdω (2.37)

2.3 soluciones fundamentales 13

2.2.2.4 Desplazamiento y tracciones

Alternativamente a la Ecuación (2.20), los desplazamiento en función de los poten-ciales pueden ser expresados en forma matricial como:

ux

uy

uz

=

∂∂x − ∂

∂z 0

0 0 ∂∂x

∂∂z

∂∂x 0

Φ

Ψ

χ

(2.38)

que muestra claramente el desacoplo de los desplazamiento en el plano y fuera delplano.

Las tracciones sobre un plano horizontal de normal ~ez pueden derivarse del cam-po de desplazamiento sin más que aplicar las Ecuaciones de compatibilidad (2.1) ycomportamiento (2.2):

tx

ty

tz

~ez

=

0 0 0 0 0 2µ

0 0 0 0 2µ 0

λ λ λ + 2µ 0 0 0

∂∂x 0 0

0 ∂∂y 0

0 0 ∂∂z

12

∂∂y

12

∂∂x 0

0 12

∂∂z

12

∂∂y

12

∂∂z 0 1

2∂

∂x

ux

uy

uz

(2.39)

2.3 soluciones fundamentales

En este capítulo se ha asentado la base matemática sobre propagación de ondas deun medio isótropo elástico y lineal en ausencia de fuerzas de volumen, y se ha dadosignificado físico a las diferentes ondas que forman el movimiento. Sin embargo elproblema particular de propagación bajo una fuerza de volumen concentrada de tipoimpulsiva tiene especial importancia desde el punto de vista ingenieril. La solución aeste problema se conoce como solución fundamental o función Green y es la base paraevaluar la respuesta dinámica de otros problemas. Por ejemplo el Método de los Ele-mentos de Contorno (MEC) es una herramienta numéricamente que permite resolverproblemas de propagación con condiciones de contorno complejas, superficies irregu-lares o incluso en medios infinitos. Éste método se basa en el teorema de reciprocidady relaciona dos estados: uno virtual conocido que corresponde con la solución funda-mental, y otro desconocido real. En los capítulos siguientes se trata en profundidadeste problema.

Matemáticamente la función de Green en desplazamiento es un tensor de segun-do orden que se escribe de la forma Gij(x, t, x′ , t′). Representa la componente i del


vector desplazamiento u en la posición x en el tiempo t debido a una carga concen-trada impulsiva en la posición x′ para el tiempo inicial t′ y en la dirección j. Estacarga se incorpora a la ecuación de Navier como una fuerza volumétrica de la forma:ρbi(x, t) = δ(x− x′)δ(t− t′)δij para la dirección ej, donde δij es la delta de Kroneckery δ e la delta de Dirac. Evidentemente la función de Green satisface la ecuación deNavier excepto en el punto x′ que existe una singularidad. La función de Green entracciones Tk

ij(x, t, x′, t′), necesario también en el MEC, tiene un significado análogo y

representa la tensión en la dirección i debido al impulso en dirección j para el planode normal unitaria nk. Cada columna se obtiene sin más que aplicar el principio deCauchy ti(x, t) = σ

jik(x, t)nk al tensor de tensiones σ

jik(x, t) derivado del campo de

desplazamiento para la carga en la dirección j.

3SOLUCIÓN FUNDAMENTAL DEL SEMIESPACIO 3DHOMOGÉNEO EN EL DOMIN IO DEL T IEMPO

3.1 introducción

En este capítulo se presenta la solución al problema de Lamb [27] presenta por John-son [19] en 1973 sobre el campo de desplazamiento de un semiespacio debido a unafuerza concentrada impulsiva. Hay numerosos artículos que abordan el mismo pro-blema [8] [9] [12] [31] [36] [37] [38] sin embargo éste trabajo destaca porque presentauna formulación realmente tridimensional (no supone invariabilidad en ninguna di-rección ni periodicidad de la fuente), y es válida para cualquier tipo de suelo y parauna fuente ubicada en cualquier posición del semiespacio así como el receptor. Ade-más proporciona expresiones explicitas del vector tracción para un plano de normalunitaria genérica nk dentro del semiespacio derivado de los desplazamientos (ambasvariables necesarias para la formulación del MEC). Sin embargo, dichas soluciones nose presentan en forma cerrada, sino en forma integral que requiere de una evaluaciónnumérica costosa por las múltiples singularidades. Éste trabajo fue ampliado poste-riormente por Triantafyllidis en su tesis doctoral [48] quien además lo utilizó comosolución fundamental del MEC.

Siguiendo la línea de trabajo de Cagniard y Hoop [11], Johnson plantea la solu-ción en el dominio transformado de Laplace, creando una dominio virtual en el quelas ecuaciones se simplifican. En este dominio se despeja el campo de desplazamien-tos para dos casos: uno particular en el que la fuente o el receptor se localiza en lasuperficie, y otro genérico en el que ambos se encuentra dentro del semiespacio. Poste-riormente la solución se revierte al dominio natural mediante un cambio de variablesingenioso conocido con el método de Cagniard-de Hoop.

En este capítulo se muestra la formulación del problema general y en el dominiotransformado. Se expone la resolución de ambos problemas, el particular y el general,se comentan el significado físico de las variables, y se explican ciertas cuestiones im-portantes de las ecuaciones y de su implementación. Concretamente se comentan losproblemas numéricos que aparecen cuando la fuente y el receptor tienden a la super-ficie, y se proponen ciertas técnicas para paliar su defecto. Por último se representanel campo de desplazamientos a lo largo del tiempo para algunos ejemplos, tras trasresolver los problemas numéricos mencionados.

15

16 solución fundamental del semiespacio 3d homogéneo en el dominio del tiempo

3.2 planteamiento del problema

Conviene en primer lugar aclarar la nomenclatura de la que se va a hacer uso a lolargo de este capítulo. Se utiliza en principio el sistema de coordenadas cartesianas, sibien, se utilizan también otras variables intermedias propias del sistema esférico. Seconsidera el semiespacio cuya superficie libre corresponde con el plano z = 0 y eje~ez normal a la superficie orientado hacia el interior del semiespacio. Se establece, sinpérdida de generalidad, que la fuente se ubica en el centro del sistema de coordena-das a la profundidad (0, 0, z′), mientras que el punto de evaluación puede ubicarselibremente en la posición genérica (x, y, z). Véase figura (3.1) para mejor aclaración.

r

b

ϴ

z

x

y

(0,0,z')

(x,y,z)

Figura 3.1: Geometría del semiespacio: superficie libre coincide con el plano z = 0 tal que el eje z positivodefine la región del semiespacio. Se evalúan los desplazamientos u en el punto (x, y, z) debido a la fuenteb ubicada en la posición (0, 0, z′)

Matemáticamente el problema en desplazamientos de un semiespacio homogéneoelástico sometido a una fuerza de volumen puede plantearse según la ecuación deNavier:

(λ + 2µ)∇∇ · u− µ∇×∇× u+ b = ρu (3.1)

siendo b el vector de fuerzas de volumen, que para el caso de una carga concentradaen la posición x′ e impulsiva en el instante t′ se expresa:

b(x, t, x′ , t′) = (bx~ex + by~ey + bz~ez)δ(x)δ(y)δ(z− z′)δ(t− t′) (3.2)

No obstante en lo que sigue se toma una fuerza dinámica de tipo escalón (H(t −t′)) que simplifica las ecuaciones. La respuesta a la carga impulsiva tipo (δ(t − t′))se calcula sin más que realizar la derivada temporal del campo de desplazamientosobtenido anteriormente. Además por comodidad, y sin pérdida de generalidad, seconsidera que la fuente actúa en el tiempo t′ = 0. De esta forma la Ecuación se redefinecomo:

b(x, t, x′) = (b1~ex + b2~ey + b3~ez)δ(x)δ(y)δ(z− z′)H(t) (3.3)

3.2 planteamiento del problema 17

con

H(t− t′) =

1 ∀t ≥ t′

0 ∀t < t′(3.4)

El problema se cierra imponiendo las condiciones iniciales y de contorno. Se suponeque el semiespacio se encuentra en reposo hasta el momento en el que se introduce laperturbación. Por otro lado, la superficie libre del semiespacio está exenta de tensiones.Matemáticamente esto se expresa como:

ui(x, 0) = 0 ∀x ∈ Ω

ui(x, 0) = 0 ∀x ∈ Ω (3.5)

ti(x, t) = 0 ∀z = 0 ∈ Γt (3.6)

La solución al problema se alcanza tras resolver el sistema de ecuaciones diferencia-les (3.1), resultando un campo de desplazamientos que de forma genérica se puedeexpresar como:

u(x, t) =

ux(x, t)

uy(x, t)

uz(x, t)

⇒ Solución al vector fuerza

bx

by

bz

δ(x)δ(y)δ(z− z′)H(t)

(3.7)Esta solución genérica ante una carga en dirección cualquiera puede expresarse

alternativamente como una combinación de soluciones básicas. Esto es, por simpleaplicación del principio de superposición puede descomponerse en la suma de las so-luciones relativas a cada una de las componentes del vector de fuerzas en la direcciónde su vector director. Dicho en forma matemática:

u(x, t) =

uxx(x, t)

uyx(x, t)

uzx(x, t)

︸︷︷︸

bx(x,t,x′)ex

+

uxy(x, t)

uyy(x, t)

uzy(x, t)

︸︷︷︸

by(x,t,x′)ey

+

uxz(x, t)

uyz(x, t)

uzz(x, t)

︸︷︷︸

bz(x,t,x′)ez

(3.8)

Y retomando el significado del tensor de Green se tiene que:

ui(x, t) = Gij(x, t, x′)bj (3.9)

Por tanto, para obtener el tensor de Green de segundo orden completo basta consustituir las componentes del vector fuerza por una matriz identidad 3x3, o lo quees lo mismo, resolver tres problemas para cada una de las direcciones cartesianas delvector fuerza con amplitud unitaria.


bj =

bx

by

bz

⇒

1 0 0

0 1 0

0 0 1

= Bj (3.10)

3.3 solución en el dominio transformado de laplace

Tal como se ha explicado, la solución en desplazamiento del problema se obtiene alresolver el sistema de ecuaciones diferenciales no homogéneo (3.1) sujeto a las condi-ciones iniciales y de contorno (3.5). De la teoría clásica sobre ecuaciones diferencialesse desprende que la solución es única, y puede plantearse como la suma de la solu-ción particular más las soluciones homogéneas del problema sin fuerzas de volumen.De un simple vistazo es evidente que este sistema de ecuaciones no es inmediato deresolver, para lo cual Johnson [19] propone transformarlo en un sistema matemática-mente más sencillo fruto de aplicar la transformada de Laplace respecto a las variables(t, x, y, z), y dando lugar a un dominio virtual homólogo (Ls,ξx ,ξy ,ξz). Luego, una vezobtenida la solución en el dominio virtual, se realiza la antitransformada al dominionatural haciendo uso del método de Cagniard and De Hoop [11].

3.3.1 Solución particular

Efectivamente, si se aplica la transformada de Laplace a la ecuación general (3.1) sereduce al siguiente sistema de ecuaciones:

ξ2x +µ

λ+µ (ξ2z − ν2s ) ξxξy ξxξz

ξxξy ξ2y +µ

λ+µ(ξ2z − ν2s ) ξyξz

ξxξz ξyξz ξ2z +µ

λ+µ (ξ2z − ν2s )

G(¸, s; 0, 0, z′ , 0) = −1

s

e−ξzz′

λ + µB

(3.11)donde G es la transformada de Laplace del tensor de Green en desplazamiento, y B

es la matriz de fuerzas volumétricas introducido en (3.10) que proviene de aplicar latransformada de Laplace a la matriz de fuerzas:

G(¸, s; 0, 0, z′ , 0) = L(ξx ,ξy,ξz,s)[G(x, t; 0, 0, z′ , 0)] (3.12)

B(¸, s; 0, 0, z′ , 0) = L(ξx ,ξy ,ξz,s)[Bδ(x)δ(y)δ(z− z′)H(t)] =1se−ξzz

′B (3.13)

En las Ecuaciones anteriores (3.11) se han tomado las siguientes definiciones:

3.3 solución en el dominio transformado de laplace 19

kp =s

cpks =

s

cs

cp =

√

λ + 2µ

ρcs =

√µ

ρ

νp = (k2p − ξ2x − ξ2y) νs = (k2s − ξ2x − ξ2y) (3.14)

Adicionalmente para asegurar la convergencia de la transformación se deben satis-facer las siguientes relaciones:

a) s ∈ ℜ y s ≥ 0

b) | ℜ[ξx ] |< kp

c) | ℜ[ξy] |< ℜ[√

k2p − ξ2x

]

d) | ℜ[ξy] |< ℜ[νp]e) ℜ[νp] ≥ 0

f ) ℜ[νs ] ≥ 0

(3.15)

Invirtiendo la transformada respecto de ξz (L−1(ξz)) y despejando el tensor de Greense llega a la solución particular en función de las variables (ξx, ξy, z, s; z′), o lo que eslo mismo a la solución en el dominio transformado Lξx ,ξy,s:

Gp(ξx, ξy, z, s; 0, 0, z′ , 0) =

ξ2x ξxξy ξxνpsign(z′ − z)

ξxξy ξ2y ξyνpsign(z′ − z)

ξxνpsign(z′ − z) ξyνpsign(z′ − z) νp)

e−νp|z′−z|

2ρs3νpB

+

ξ2x + ν2s −ξxξy −ξxνssign(z′ − z)

−ξxξy ξ2y + ν2s −ξyνssign(z′ − z)

−ξxνssign(z′ − z) −ξyνssign(z′ − z) ξ2x + ξ2y

e−νs|z′−z|2ρs3νs

B (3.16)

3.3.2 Solución homogénea

Por otro lado, la solución homogénea se calcula tomando el vector de fuerzas nuloen la Ecuación de Navier (3.1). Si se analiza detenidamente, este caso correspondecon el problema genérico de propagación de ondas en un semiespacio en ausencia


de fuerzas estudiado en el Capítulo 2. Su resolución puede abordarse mediante elmétodo de los potenciales, en el que el vector desplazamiento se descompone en unaparte irrotacional y otra solenoidal:

u = ∇Φ +∇× Ψ (3.17)

Introduciéndolo en la ecuación de Navier se llega nuevamente a un sistema de dosecuaciones desacopladas en función de potencial φ y del vector Ψ.

∇2Φ =1c2p

Φ (3.18)

∇2Ψ =

1c2s

Ψ (3.19)

Ahora, siguiendo la propuesta de Johnson [19] y tomando la transformada de Laplacerespecto de las variables (x, y, t), se llega a dos ecuaciones diferenciales ordinarias desegundo orden según la coordenadas vertical z:

∂2Φ

∂z2− ν2pΦ = 0 (3.20)

∂2Ψ

∂z2− ν2s Ψ = 0 con Φ, Ψ = Lξx ,ξy,s[Φ,Ψ] (3.21)

En el que teniendo en cuenta que el semiespacio está definido en la region positiva deleje z, la solución de sendos potenciales se define:

Φ = φ(ξx, ξy, s)e−νpz) (3.22)

Ψ = Œ(ξx, ξy, s)e−νpz (3.23)

Por tanto, recomponiendo el vector desplazamiento, la solución homogénea en eldominio transformado Lξx ,ξy,s es:

Gh = G1h + G2

h =

a1ξx a2ξx a3ξx

a1ξy a2ξy a3ξy

−a1νp −a2νp −a3νp

e−νp(z′−z)

2ρs3νpB

b1νs b2νs b3νs

c1νs c2νs c3νs

b1ξx + c1ξy b2ξx + c2ξy b3ξx + c3ξy

e−νs(z′−z)

2ρs3νsB (3.24)

Donde ai, bi, ci son coeficientes que se deducen de las condiciones de iniciales y decontorno, y en el que ya se ha tenido en cuenta la condición ∇ · Ψ = 0.

3.3 solución en el dominio transformado de laplace 21

3.3.3 Solución en el dominio transformado

Obviamente, uniendo ahora la solución particular (3.16) y la solución homogénea(3.24) se alcanza finalmente la solución completa en el dominio transformado Lξx ,ξy ,s

. Si se observa, esta solución está formada por cuatro términos, dos relativos a lasolución particular y otros dos términos de la solución homogénea, carentes en unprincipio de significado físico. Reordenando, se puede expresar alternativamente:

G(ξx, ξy, z, s; 0, 0, z′) =[P(ξx, ξy, z, s; 0, 0, z′) + S(ξx, ξy, z, s; 0, 0, z′) + PP(ξx, ξy, z, s; 0, 0, z′)

+SS(ξx, ξy, z, s; 0, 0, z′) + PS(ξx, ξy, z, s; 0, 0, z′) + SP(ξx, ξy, z, s; 0, 0, z′)]B (3.25)

siendo:

P(ξx, ξy, z, s; z′) =



ξxνpsign(z′ − z) ξyνpsign(z′ − z) ν2p

e−νp|z′−z|

2ρs3νp

(3.26a)

S(ξx, ξy, z, s; z′) =

ξ2y + ν2s −ξxξy −ξxνssign(z′ − z)

−ξxξy ξ2y + ν2s −ξyνssign(z′ − z)

−ξxνssign(z′ − z) −ξyνssign(z′ − z) ξ2x + ξ2y

e−νs|z′−z|

2ρs3νs

(3.26b)

PP(ξx, ξy, z, s; z′) =



ξxνp ξyνp ν2p

be−νp(z′−z)

2ρs3νpd(3.26c)

SS(ξx, ξy, z, s; z′) =

(ξ2x + ν2s )d− 8νpν3s ξ2x −ξxξy(d+ 8νpν3s ) −ξxνsb

−ξxξy(d+ 8νpν3s ) (ξ2x + ν2s )d− 8νpν3s ξ2y −ξyνsb

ξxνsb ξyνsb −(ξ2x + ξ2y)b

e−νs(z′−z)

2ρs3dνs

(3.26d)


PS(ξx, ξy, z, s; z′) =

ξ2xνs ξxξyνs ξxνpνs

ξxξyνs ξ2yνs ξyνsνp

ξx(ξ2x + ξ2y) ξy(ξ2x + ξ2y) νp(ξ2x + ξ2y)

2h

e−(νpz′+νsz)

ρs3d(3.26e)

SP(ξx, ξy, z, s; z′) =

ξ2xνs ξxξyνs −ξx(ξ2x + ξ2y)

ξxξyνs ξ2yνs −ξy(ξ2x + ξ2y)

−ξxνpνs −ξyνsνp νp(ξ2x + ξ2y)

2h

e−(νpz′+νsz)

ρs3d(3.26f)

con

h = nu2s − ξ2x − ξ2y

d = h2 + 4νsνp(ξ2x + ξ2y)

b = h2 − 4νsνp(ξ2x + ξ2y) (3.27)

que corresponde con la solución presentado por Johnson [19]. Tal como sugiere lanotación, ahora cada término sí tiene un sentido físico, y representa el movimientode las ondas P que alcanzan el receptor directamente sin reflexión alguna (P), de lasondas transversales directas (S), de las ondas P reflectadas en la superficie (PP), delas ondas S reflectadas (SS), de las ondas P reflectadas que cambian su naturaleza aondas S (PS) y viceversa(SP), respectivamente.

3.4 respuesta sobre la superficie libre

Con la idea de avanzar ordenadamente desde el caso más sencillo hasta el más ge-nérico, en este apartado se muestra la solución al problema particular en el que elreceptor se encuentra sobre la superficie libre (z = 0), mientras que la fuente puedeubicarse en cualquier posición dentro dentro del semiespacio. Es evidente que éstasuposición simplifica el problema y promueve la reagrupación de las matrices funda-mentales. Cabe destacar aquí que, si bien las matrices que se manejan son más sencillasque para el caso general, debe entenderse que el modo de proceder es totalmente váli-do para éste último.Efectivamente, tomando z = 0 en la solución general en el dominio transformado

(Ec. (3.25)) la expresión se simplifica, de modo que puede escribirse únicamente enfunción de dos términos tal como se muestra a continuación:

3.4 respuesta sobre la superficie libre 23

G(ξx, ξy, 0, s; 0, 0, z′) =

M(ξx, ξy, 0, s; 0, 0, z′)e−νpz

′

µsdB+ N(ξx, ξy, 0, s; 0, 0, z′)

e−νsz

µsdB (3.28)

Ahora físicamente el primer término aglomera la parte dilatacional del movimientoy el otro recoge el comportamiento de las ondas transversales.

El siguiente paso consiste en aplicar la transformada inversa de Laplace al dominionatural L−1(ξx ,ξy,s), pero evidentemente no es inmediato. En este punto, primero Cag-niard y despues de Hoop [11] desarrollaron un ingenioso método que por medio deuna serie de sustituciones de variables que hacen posible la transformación L−1(ξx ,ξy).La transformación al tiempo puede realizarse por simple inspección, tras aplicar elmétodo, como se demostrará. A continuación se muestra la esencia del método que esde aplicación a cualquier término de la Ecuación (3.28) M, N o de la Ecuación general(3.25. De hecho se toma un término genérico que se notará por G:

1. Se pretende realizar la transformada inversa L−1(ξx ,ξy)[G]. Como paso previo sesepara del término G la parte exponencial en función de la profundidad (z, z′) yde la variable ν:

G(x, y, z, s; 0, 0, z′) =−14π2

i∞x

−i∞

[

G∗(ξx, ξy, s)e−νγ f (z,z′)]

eξxx+ξyydξxdξy (3.29)

donde:

νγ = νp o νs según el término se refiera a las ondas de compresión o cizalladura

G∗ ≡ cualquier término M, N , o de la solución general P, S, PP, SS, PS, SP

sin la parte exponencial e−νγ f (z,z′) (3.30)

2. Cambio las variables de integración ξx, ξy por las variables p, q según la regla:

ξx = sq cos ϕ − isp sin ϕ

ξy = sq sin ϕ + isp cos ϕ

(3.31)

Por tanto realizando el Jacobiano de la trasformación se tiene la relación:

dξxdξy = is2dpdq (3.32)


Y desarrollando la parte exponencial en función de p y q:

e−νγ f (z,z′)+ξxx+ξyy = e−s(ηγ f (z,z′)−qR) (3.33)

donde se define:

ηγ =

ηp =√

c−2p − q2 + p2 ℜ[ηp] ≥ 0

ηs =√

c−2s − q2 + p2 ℜ[ηs ] ≥ 0

(3.34)

y se ha hecho uso de la relación geométrica:

R =√

x2 + y2 (3.35)

3. Observar que, tras el cambio de variables, la integración respecto a la variablep se define solo a la largo de su eje real. Además, teniendo en cuenta que elintegrando en función de p es simétrico respecto al eje positivo/negativo, laintegral puede reducirse a:

G(x, y, z, s; 0, 0, z′) =−24π2

∫ ∞

0

∫ i∞

−i∞G∗(p, q)e−s(ηγ f (z,z′)−qR)dqdp (3.36)

4. Se introduce la variable τγ que se entiende como un tiempo adimensional:

τγ = ηγ f (z, z′)− qR τγ =

τp

τs(3.37)

Por otro lado, notar que la función f (z, z′) geométricamente expresa la distanciavertical recorrida por las ondas a las que el término en cuestión hace referencia,y que por supuesto incluyen la reflexión en su caso. De este modo, si se defineun nuevo radio vector que describe la distancia total recorrida, a modo similarque el método de las imágenes, se tiene que:

f (z, z′) = r cos θ

R = r sin θ (3.38)


Sustituyendo en la definición de τ:

τγ = ηγ r cos θ − qr sin θ (3.39)

5. La definición de la variable temporal τγ permite estudiar el comportamiento dela variable de integración q en función de ésta, deduciéndose que en el planoq-complejo tiene forma de hipérbola simétrica respecto del eje real.

τγ =√

c−2γ − q2 + p2r cos θ − qr sin θ

(3.40)

Despejando simplemente:

qγ = −τγ

rsin θ ± ir cos θ

√(τγ

r

)2− c−2

γ − p2 (3.41)

que para un valor dado de la variable de integración p, se representa en el planocomplejo en función de τγ tal como se muestra en la figura (3.2)

Imq

Req

Vértice de la hipérbola:

A B

CD

Figura 3.2: Representación de la variable de integración q en el plano complejo. Esta representación tomael nombre de Contorno de Cagniard. Triantafyllidis [48]


En lo que se refiere al camino de integración de la variable q (−i∞, i∞), aplicandoel teorema integral de Cauchy sobre un contorno cerrado se cumple la relación:

∫ C

BI +

∫ D

CI +

∫ A

DI +

∫ B

AI = 0 (3.42)

Siendo I el integrando en función de la variable q de la ecuación (3.36). Y teniendoen cuenta que:

∫ D

CI +

∫ B

AI = 0 (3.43)

Se deduce que la integral a lo largo del eje imaginario es equivalente a la integrala lo largo de la curva qγ:

∫ C

BI =

∫ D

AI = 0

⇒∫ i∞

−i∞I =

∮

I|qγ (3.44)

6. Se sustituye la variable q por la variable temporal τ. Notar además que el inte-grando es simétrico respecto al eje real de q, luego solo permanece en la integralla parte imaginaria.

Es fácil deducir la relación:

dqγ =

− sin θ

r± i

1√( τγ

r

)2 − c−2γ − p2

τγ

rcos θ

dτγ (3.45)

Y para el caso particular del vértice de la hipérbola:

τ∗γ = r

√

c−2γ + p2 (3.46)

Por tanto, desarrollando la integral (3.36) y teniendo en cuenta todos los pasosse llega a:


G(x, y, z, s; 0, 0, z′) =

−24π2

∫ ∞

0

∮ ∞

∞G∗(p, q, τγ)|qγ

− sin θ

r± i

1√( τγ

r

)2 − c−2γ − p2

τγ

rcos θ

e−sτγdτγdp

=1

π2

∫ ∞

0ℜ

∫ ∞

r√

c−2γ +p2

G∗(p, q, τγ)|q+γ

1

√( τγ

r

)2 − c−2γ − p2

τγ

rcos θ

e−sτγ

dτγdp

(3.47)

que de forma compacta puede expresarse:

G(x, y, z, s; 0, 0, z′) =1

π2

∫ ∞

0ℜ∫ ∞

r√

c−2γ +p2

F(p, q, τγ, ...)e−sτγ

dτγdp (3.48)

7. Por último se realiza la transformada inversa temporal L−1(s). De la Ecuaciónanterior en forma compacta (3.48) se deduce que la anti-transformada de Laplacees inmediata:

G(x, y, z, t; 0, 0, z′) = L−1(s)G(x, y, z, s; 0, 0, z′)

=1

π2

∫ ∞

0ℜ

G∗(p, q, τγ)|q+γ

1

√( τγ

r

)2 − c−2γ − p2

τγ

rcos θ

H

(

t− r√

c−2γ + p2

)

dp

(3.49)

y simplificando:

G(x, y, z, t; 0, 0, z′) =

=1

π2

∫√

(t/r)2−c−2γ

0ℜ

G∗(p, q, t)|q+γ

ηγ

√(tr

)2 − c−2γ − p2

H(

t− rc−1γ

)

dp

(3.50)

Ahora bien, si se particulariza el método de la transformada inversa de Cagniard-deHoop para el problema particular (Ec. (3.28)), se alcanza la solución en el dominioespacial-temporal presentado por Johnson [19]. El resultado final tras sustituir los tér-minos es tal como se expresa:


G(x, y, 0, t; 0, 0, z′) =

=1

π2µr

∫√

(t/r)2−c−2p

0ℜ

ηp

√(tr

)2 − c−2p − p2

M(p, q, t)

σ

H(

t− rc−1p

)

dpB

=1

π2µr

∫ p2

0ℜ

ηs

√(tr

)2 − c−2s − p2

N(p, q, t)

σ

H (t− t2) dpB (3.51)

Siendo:

M11 = 2ηs

((q2 + p2) cos2 ϕ − p2

)

M12 = 2ηs(q2 + p2) sin ϕ cos ϕ

M13 = 2qηpηp cos ϕ

M21 = 2ηs(q2 + p2) sin ϕ cos ϕ

M22 = 2ηs

((q2 + p2) sin2 ϕ − p2

)

M23 = 2qηpηp sin ϕ

M31 = qγ cos ϕ

M32 = qγ sin ϕ

M33 = ηpγ

(3.52)

N11 = η−1s

[η2s γ − (γ − 4ηsηp)((q2 + p2) sin2 ϕ − p2)

]

N12 = η−1s (γ − 4ηsηp)(q2 + p2) sin ϕ cos ϕ

N13 = −qγ cos ϕ

N21 = η−1s (γ − 4ηsηp)(q2 + p2) sin ϕ cos ϕ

N22 = η−1s

[η2s γ − (γ − 4ηsηp)((q2 + p2) cos2 ϕ − p2)

]

N23 = γ sin ϕ

N31 = −2qηpηs cos ϕ

N32 = −2qηpηs sin ϕ

N33 = 2ηp(q2 − p2)

(3.53)

γ = η2s + p2 − q2 (3.54)

σ = γ2 + 4 ∗ ηpηs(q2 − p2) (3.55)


p2 =

√

(t/r)2 − c−2s sin θ ≤ cs/cp

√(

t/r−√

c−2s −c−2

p cos θ

sin θ

)2

− c−2p sin θ > cs/cp

(3.56)

t2 =

r/cs sin θ ≤ cs/cpr/cp sin θ + r

√

c−2s − c−2

p cos θ sin θ > cs/cp(3.57)

y:

q =

− tr sin θ + ir cos θ

√(tr

)2 − c−2p − p2 para el término M

− tr sin θ + ir cos θ

√(tr

)2 − c−2s − p2 para el término N

(3.58)

Figura 3.3: Representación del semiespacio para el problema particular en el que receptor se localizasobre la superficie

Se debe notar que aquí se ha expresado la solución al problema en el que la fuen-te se encuentra dentro del semiespacio y el punto de evaluación sobre la superficie.Sin embargo esta solución es equivalente al problema inverso, en el que la fuente selocaliza sobre la superficie y el receptor a una profundidad z′. Sin más que aplicar larelación de reciprocidad se demuestra:

Gij(x, y, 0, t; 0, 0, z′ , t′) = Gji(−x,−y, z′, t− t′; 0, 0, 0, 0) (3.59)


3.4.1 Resultados numéricos

La función de Green que se ha presentado en esta sección (caso particular del re-ceptor sobre la superficie) ha sido programada en Matlab para ser evaluada numérica-mente, y con la idea que pueda servir como solución fundamental del MEC. En estasección se comentan los problemas numéricos que se han encontrado y se explicanformas de proceder para evitarlos o al menos paliarlos. Además se han procesado cier-tos problemas sencillos y se han validado los resultados mediante comparación con lasolución explícita de Pekeris [38].Tal como se ha presentado, todos los términos de la función de Green están ex-

presados en forma integral, lo cual no es un problema desde el punto de vista de laimplementación puesto que existen numerosas técnicas para evaluarlas, pero el nú-cleo contiene numerosas singularidades. Si se estudia el término correspondiente alas ondas de compresión M, es fácil percatar que tiene una singularidad cuando lavariable de integración se acerca al límite de integración superior ((t/r)2 − 1/c2p)

1/2.Dicha singularidad es débil, esto es, aunque el integrando tiende a infinito a medidaque se aproxima al polo, el área que encierra la curva es finita, y por tanto es integra-ble. Dicho matemáticamente, el integrando tiene la forma de 1/(x− x0)n en el que elexponente es menor que la unidad. De hecho, Matlab incluye un paquete de funcionespara realizar integrales numéricas, entre ellas destacar la función quadgk que se basaen la cuadratura Gauss-Kronrod adaptativa [45] y que es capaz de hallar integralessingulares débiles. Aún así, es conveniente evitar estas singularidades en la medida delo posible puesto que conlleva un alto coste computacional. En este caso Triantafyllidis[48] propuso un simple cambio de variables tal como se demuestra a continuación.

M(x, t, x′) =∫√

(t/r)2−1/c2p

0

F(p, q, z, z′)

√

(t/r)2 − 1/c2p − p2

dp

=∫√

(t/r)2−1/c2p

0

F(p, q, z, z′)√√√√

p−√

(t/r)2 − 1/c2p︸︷︷︸

p1

√√√√

p+√

(t/r)2 − 1/c2p︸︷︷︸

p2

dp

∼∫ p1

0

1√p− p1

dp =⇒ singularidad débil (3.60)

Cambio de variable:

p =√

(t/r)2 − 1/c2p sin v (3.61)


Sustituyendo en la integral puede observarse que la singularidad desaparece:

M(x, y, 0, t; 0, 0, z′) =1

π2µr

∫ π/2

0ℜ(ηp

σ

)

M(p, q, t)

H(

t− rc−1p

)

dv

(3.62)

Por otro lado, al término correspondiente a las ondas transversales N le sucede prác-ticamente lo mismo: presenta una singularidad débil. Pero en este caso la singularidadpuede hallarse justamente en el límite superior de integración o en medio del caminode integración, depende de la relación entre las velocidades de las ondas S y P. Sedistinguen dos casos:

El tiempo de llega de las ondas S es menor que el de las ondas SP (sin θ < cs/cp):En este caso la singularidad se localiza en el extremo superior del límite deintegración, por lo que se procede igual que para el término M, es decir, se

realiza el cambio de variable p =√

(t/r)2 − 1/c2p sin v quedando un resultado

análogo.

El tiempo de llegada de las ondas SP es menor que el de las ondas S. En estecaso aparece una singularidad en medio del camino de integración sólo si eltiempo de evaluación es mayor que el tiempo de llegada de las ondas S, es decir,t > r/cs. Entonces la integral debe separarse en dos y aislar la singularidad enlos extremos. Para cada integral Triantafyllidis [48] propuso un cambio diferentede variables:

N(x, y, 0, t; 0, 0, z′) =1

π2µr

∫ p2

0ℜ

ηs

√(tr

)2 − c−2s − p2

N(p, q, t)

σ

H (t− t2) dpB

=1

π2µr

∫√

(t/r)2−1/c2s

0F(p, q, t, z, z′)dp

︸︷︷︸

I1

+1

π2µr

∫ p2√

(t/r)2−1/c2sF(p, q, t, z, z′)dp

︸︷︷︸

I2

(3.63)

Para la integral I1 se realiza el mismo cambio de variable que se ha utilizado

para el caso anterior p =√

(t/r)2 − 1/c2p sin v, mientras que para la integral I2

se realiza un cambio alternativo p = v2 +√

(t/r)2 − 1/c2s . El resultado es:


N(x, y, 0, t; 0, 0, z′) =1

π2µr

∫ π/2

0ℜηs

σN(p, q, t)

H (t− t2) dvB

+1

π2µr

∫ 0√

p2−√

(t/r)2−1/cs2ℜ

ηs

σ

2√

−v2 − 2√

(t/r)2 − 1/cs2N(p, q, t)

H (t− t2) dvB

(3.64)

Después de las transformaciones comentadas, las integrales quedan exentas de sin-gularidades. En las Figuras 3.4-3.7 se muestran algunos ejemplos numéricos de laevolución temporal de los desplazamientos de un semiespacio rocoso con las mis-mas propiedades que las elegidas previamente por Johnson [19] y Triantafyllidis [48](cp = 8000 m/s, cs = 4620 m/s, ρ = 3300 kg/m3). Concretamente se muestran lostérminos no nulos del tensor de Green para tres configuraciones diferentes: la primeracorresponde al caso de un baja ratio de la distancia horizontal entre la profundidad(Figura 3.4); La siguiente configuración se refiere a una situación inversa (Figura 3.5);Mientras que la tercera configuración intenta representar un caso límite en el que lafuente tiende a la superficie (Figura 3.6 y Figura 3.7). Para evaluar las integrales seha utilizado la función quadgk que tiene implementada Matlab, y que utiliza el mé-todo de la cuadratura adaptativa Gauss-Kronrod. Se han evaluado los desplazamien-tos hasta un tiempo final igual a t f inal = 2r/cs y se ha escogido un paso temporal∆t = t f inal/500. Las gráficas se muestran en forma adimensional: el eje de ordenadas,que representa cada una de las componentes de la función de Green adimensionali-zadas por GH

ij = Gijµπr/Z siendo Z la magnitud de la fuerza; mientras que el eje deabscisas temporal por τ = tcs/rEn primer lugar se puede observar que culitativamente cualitativamente que los

resultados obtenidos están en consonancia con los mostrados por Johnson [19] Porotro lado, de los resultados pueden extraerse varias observaciones:

En todos los ejemplos mostrados se ha localizado la fuente y el receptor en elplano y = 0. Por tanto los desplazamientos en el plano (ux, uz) están originadospor las ondas P-SV, mientras que los desplazamientos fuera de él (uy) se debena las ondas horizontalmente polarizadas SH.

Se aprecia que la respuesta depende en gran medida del ángulo de emisión θ, olo que es lo mismo, de la relación entre la distancia horizontal y la profundidadrelativa fuente-receptor. Para ángulos de emisión pequeños, que corresponde conla Figura 3.4, las ondas directas P y S son las principales responsables del cam-po de desplazamientos. En el caso de ángulos de emisión altos, correspondientecon la Figura 3.5, aparecen nuevos tipos de ondas que son las reflejas SP y lasde Rayleigh. Éste último tipo se propagan por la superficie y se atenúan rápi-damente con la profundidad. En la Figura 3.6 se muestra un ejemplo en el que


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

P S

(a) GHxx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) GHxz

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) GHyy

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) GHzx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(e) GHzz

Figura 3.4: Evolución temporal de las componentes de la función Green adimensional, GH(2 ×103 , 0, 0, t; 0, 0, 10× 103), de un semiespacio homogéneo. Ángulo de emisión θ = 11.31o.

la fuente y el receptor están próximos a la superficie y las ondas de Rayleighson las dominantes. Y en la Figura 3.7 se muestra un caso límite en el que lafuente prácticamente toca la superficie y las ondas de Rayleigh produce que losdesplazamientos tiendan a infinito, esto es, aparece una singularidad.


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

P SP SR

(a) GHxx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) GHxz

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) GHyy

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) GHzx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(e) GHzz

Figura 3.5: Evolución temporal de las componentes de la función Green adimensional, GH(10 ×103, 0, 0, t; 0, 0, 2× 103), de un semiespacio homogéneo. Ángulo de emisión θ = 78.69o.

Se encuentra inestabilidad numérica a partir del tiempo de llegada de las ondasde Rayleigh cuando la fuente está muy próxima a la superficie (3.7). Ello se debea que se excitan los polos de Rayleigh reales (polos de la ecuación σ(p) = 0) y losintegrandos se vuelven singulares. Cuando el ángulo de emisión tiende a pi/2


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

P SR

(a) GHxx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) GHxz

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) GHyy

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) GHzx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(e) GHzz

Figura 3.6: Evolución temporal de las componentes de la función Green adimensional, GH(10 ×103 , 0, 0, t; 0, 0, 0.2× 103), de un semiespacio homogéneo. Ángulo de emisión θ = 88.85o.

los polos de Rayleigh complejos se desplazan hacia el eje p-real y afectan a laintegral. Pero evidentemente la solución debe existir, por lo que la integral tieneque tener solución finita desde el punto de vista del Valor Principal de Cauchy


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

P SR

(a) GHxx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) GHxz

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) GHyy

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) GHzx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(e) GHzz

Figura 3.7: Evolución temporal de las componentes de la función Green adimensional, GH(10 ×103, 0, 0, t; 0, 0, 1× 10−5), de un semiespacio homogéneo. Ángulo de emisión θ ≃ 90o.

(VPC). A continuación se explican algunas técnicas numéricas para paliar esteproblema.

En la Figura 3.8 se muestra el tipo de singularidad que posee el integrando en el ejereal de p cuando la fuente y receptor se encuentra sobre la superficie. Concretamente


se muestra la singularidad del integrando para el término Mxx, pero es análoga parael resto de término. Es una singularidad fuerte no integrable. Está claro que la integraldebe tener solución desde el puno de vista del VPC, y en la Referencia [34] se detallantécnicas para alcanzar la solución. Entre ellas destacar el método muy extendido dela regularización, que consiste en dividir la integral en dos, de forma que una deellas se reduce a una regular (puede evaluarse por simple cuadratura) y la restantepuede ser resulta analíticamente. O la técnica del cambio de variable que elimina lasingularidad directamente. En este sentido, se debe notar que el núcleo de la integraltiene una expresión muy compleja, incluso cambia para cada término de la funciónde Green, y no tiene forma de fracciones polinómicas. Por tanto cualquiera de losmétodos analíticos o semi-analíticos comentados son de difícil aplicación. Entonces, seentiende que debe abordarse numéricamente.

Figura 3.8: Tipo de singularidad del integrando del término M para ángulo de emisión θ = π/2

En primer lugar es lógico estudiar los polos de Rayleigh que son los responsablesde la singularidad. La expresión del término en función de p es como se expresa acontinuación:

σ(p, q) = (1/c2s − 2(q2 − p2))2 + 4(q2 − p2)√

1/cs2 − (q2 − p2)√

1/c2p − (q2 − p2) = 0

q(p) = −t sin θ/r+ i cos θ√

(t/r)2 − 1/c2γ − p2 (3.65)

Sin embargo, puede simplificarse agrupando términos y renombrando, llegándosea la función de Rayleigh:


R(ξ2) = (1− 2ξ2)2 + 4ξ2√

1− ξ2√

a2 − ξ2 = 0 (3.66)

donde se define:

ξ2 = (q2 − p2)c2s ≡ c2s/c2 → siendo c la velocidad de las ondas

a =cscp

=1− 2ν

2− 2ν(3.67)

Que operando puede transformarse en la ecuación bicúbica siguiente:

1− 8ξ2 + 8(3− 2a2)ξ4 + 16(−1+ a2)ξ6 = 0 (3.68)

La solución son tres raíces [ξ21, ξ22, ξ

33], los dos primeros no tienen significado físico,

mientras que la tercera ξ3 = cs/cR es la verdadera raíz de la ecuación. Si se estudiadetenidamente existe un valor de transición ν0 = 0.2631 por debajo del cual, las tresraíces son reales y satisfacen la relación 0 < ξ21 < ξ22 < a2 < 1 < ξ23. Si ν es exactamenteigual a ν0 entonces las raíces falsas son coincidentes ξ1 = ξ2, y si es mayor entonceslas raíces falsas son complejas y conjugadas, mientras la raíz verdadera ξ3 es siemprereal. Deshaciendo ahora el cambio de variable ξ2 = (q2 − p2)c2s se obtiene finalmenteel valor de la variable de integración p que anula la función de Rayleigh σ:

p2R =cot2 θ

c2γ+ (−ξ2/c2s + t2/r2) csc2 θ ±

2t cos θ csc2 θi√

c2γξ2/c2s − 1

cγr(3.69)

Se puede observar que solamente para θ = π/2 el valor de pR que anula la funciónde Rayleigh es real y doble, para otro caso existen dos soluciones complejas conjugadas.Además, depende del tiempo de evaluación y solamente para tiempos mayores que elde llegada de las ondas de Rayleigh (t > r/cR) el valor es positivo y dentro del caminode integración.A la luz del estudio desarrollado, y a sabiendas que los métodos semi-analíticos no

pueden aplicarse aquí, se proponen dos alternativas de paliar el efecto de los polos dela función de Rayleigh:

El primero se basa en la utilización adecuada de la potencia de la función deintegración quadgk. Como se ha comentado antes, y se recalca de nuevo, ésta escapaz de hallar integrales con singularidades débiles, es decir, las que tienen unvalor finito, si y solo si la singularidad se localiza en el extremo de la integral.La cuestión consiste en dividir la integral en dos y aislar justamente la singulari-dad en los extremos. Entonces la integral converge siempre y cuando no se tome


un valor del ángulo de emisión exactamente igual a π/2, sino suficientementepróximo para que prácticamente se acepte que la fuente es superficial. Es decir,numéricamente el nucleo no es puramente singular, y al llevar el punto proble-mático a los extremos del integrando Matlab interpreta que es una singularidaddébil que puede resolver adaptativamente. Este método, como puede verse espuramente numérico, y tiene el inconveniente que tiene un coste computacionalelevado.

El segundo método se basa en la aplicación del Teorema de los Residuos para ha-llar el VPC de la integral. La idea es interpretar el integrando como una funciónde variable compleja, e integrar esta función sobre una curva de Jordan adecuadaaplicando el teorema de los residuos. La cuestión ahora es elegir una curva cerra-da adecuada de forma que contenga el camino de integración que se pretendeaveriguar, y la integral sobre el camino restante sea sencilla de calcular. En estesentido, se propone tomar un camino complejo tal como se muestra en la Figura3.9, y de forma que rodee al polo de Rayleigh. Entonces la integral se reduce auna integral de camino complejo cuyo integrando está exento de singularidades,más la suma del residuo evaluado en el polo. Esta nueva integral puede ser eva-luada numéricamente mediante la función quadgk que permite introducirle uncamino complejo. Matemáticamente ello se expresa:

Imp

Rep

pR

p'R

p2

arco

Figura 3.9: Curva de Jordan en el plano complejo p

VPC∫ p20 F(p)dp = −

∮

arco F(p)dp+ πiRes [F(p), pR] Im[pR] = 0

VPC∫ p20 F(p)dp = −

∮

arco F(p)dp+ 2πiRes [F(p), pR] Im[pR] > 0(3.70)


Donde F(p) es el integrando de cualquiera de los términos M(x, y, 0, t; 0, 0, z′) oN(x, y, 0, t; 0, 0, z′).

Notar que mediante este método se elimina la singularidad del integrando cuan-do el ángulo de emisión es π/2, e incluso se evita pasar cerca de la singularidadpara ángulos de emisión próximos al crítico. Es decir, al desviar el camino de inte-gración evitando la singularidad o su influencia, la integral se mantiene siempreestable a un bajo coste computacional.

En la Figura 3.10 se muestra una comparativa de ambos métodos aquí explicados.Se ha tomado las mismas propiedades dinámicas que las elegidas previamente, yse ha representado la evolución temporal de los desplazamientos en la direcciónvertical debido a un carga escalón vertical. El método 1 se refiere al procedimien-to de aislar la singularidad en los extremos, mientras que le método 2 se basaen el Teorema de los residuos. En este último caso se ha tomado un ángulo deemisión exactamente igual a pi/2. El tiempo de computo hallado para el método2 ha sido el quince por ciento del método 1 y muestra una mejor estabilidadnumérica.

Ahora, una vez resulto todos los problemas numéricos, se está en condiciones devalidar la solución de Johnson por comparación con la solución explícita de Pekeris[38]. Caber recordar que Pekeris resuelve los desplazamientos superficiales en el planodebido a una carga superficial normal tipo escalón solamente para el caso de un semi-espacio homogéneos con coeficiente de Poisson de 0.25. En las Figuras 3.11 a) y 3.11 b)se muestran la evolución temporal de los desplazamientos en el plano ux(t), uz(t) res-pectivamente, calculados mediante la formulación de Johnson y Pekeris. Se ha elegidoun semiespacio con las propiedades: cs = 250 m/s, cp = 433.0127 m/s, ρ = 1850 kg/m3

Efectivamente, puede apreciarse que producen el mismo campo de desplazamientos yreproducen el efecto de las ondas S, P y de Rayleigh.Como última cuestión de este apartado, recordar que se ha obtenido la respuesta

del semiespacio homogéneo ante una carga tiempo Heaviside. La respuesta al impulsotipo delta de Dirac se alcanza por simple derivación temporal, esto es:

Gδij(x, t, x

′) =∂GH

ij (x, t, x′)

∂t(3.71)

Numéricamente la deriva temporal puede calcularse por simples diferencias, sin em-bargo este procedimiento suele tener problemas de estabilidad y precisión. Una formade evitarlo es por simple aplicación de la segunda integral de Duhamel, mediante lacual se puede obtener la solución a cualquier tipo de excitación sin más que realizaruna convolución temporal (⊗). Matemáticamente se expresa así:


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

Metodo 1Metodo 2

Figura 3.10: Comparativa de la función de Green GH33(10 × 103, 0, 0, t; 0, 0, 0) computados mediante el

método del aislamiento de la singularidad en los extremos de la integral (Método 1), y el método basadoen la aplicación del Teorema de los residuos (Método 2).

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

PekerisJohnson

(a) GH13(10× 1r3, 0, 0, t; 0, 0, o)

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

PekerisJohnson

(b) GH33(10× 1r3, 0, 0, t; 0, 0, o)

Figura 3.11: Validación de la solución de Johnson con la solución superficial en el plano de Pekeris. a)Evolución temporal del desplazamiento superficial horizontal. b) Evolución temporal del desplazamientosuperficial vertical

Uij(x, t, x′) =∂GH

ij (x, t, x′)

∂t⊗ E(t) = GH

ij (x, t, x′)⊗ ∂E(t)

∂t=

GHij (x, t, x

′)E(0) +∫ t

0GHij (x, t− τ, x′)

∂E(τ)

∂τdτ (3.72)

Siendo E una excitación genérica y Uij es el campo de desplazamiento fruto de laexcitación.

Con el objetivo de evaluar este procedimiento, se ha calculado la respuesta ante unaexcitación sencilla tipo coseno impulsivo de periodo tp. Observar que para el caso enel que el periodo de la excitación tiende a cero se reproduce la excitación impulsiva


delta de Dirac pero con un desfase igual a tp/2. En la Figuras 3.12 se muestra laevolución temporal adimensional de los desplazamientos superficiales en direcciónhorizontal de un semiespacio debido a una carga impulsiva superficial en la mismadirección (Gδ

xx). Se ha calculado de tres formas diferentes: Por simple diferencia de lasolución de Green ante la carga escalón; Mediante segunda integral de Duhamel parauna fuente tipo coseno impulsivo de duración tp = r/(120cs) computada de formadiscreta; y de la misma manera que la precedente pero computada en forma continuamediante integral numérica. Nótese, que el periodo es suficientemente pequeño talque prácticamente es un excitación tipo delta de Dirac. Las propiedades dinámicas delsemiespacio son: cs = 250 m/s, cp = 450 m/s, ρ = 1850 kg/m3. Tal como se muestra enla figura parece ser que las tres respuestas se superponen prácticamente, sin embargoobservar que la respuesta discreta de Duhamel no alcanza un permanente igual a cero,sino existe un pequeño desfase. En términos de precisión se puede concluirse que lasolución más precisa es la computada por las diferencias, seguida por la respuesta deDuhamel continua, que es tanto más precisa cuanto menor sea el periodo tp, y porúltimo las respuesta de Duhamel discreta. Sin embargo ésta última se muestra la másestable puesto que no precisa de integración numérica ni de diferencias. En lo quese refiere al coste computacional, el método de Duhamel continuo es excesivamentecostoso, y con diferencia el método de las diferencias es el más rápido.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2−250

−200

−150

−100

−50

0

50

100

150

200

250

Tiempo [−]

Des

plaz

amie

nto

[−]

DiferenciasDuhamel discretoDuhamel continuo

Figura 3.12: Evolución temporal de la función de Green Gδ11(15, 0, 0, t; 0, 0, 0) computados por tres formas:

Mediante diferencias; mediante integral de Duhamel para una excitación coseno impulsiva discreta; ypor el mismo procedimiento continuo

3.5 respuesta dentro del semiespacio

Continuando con la resolución del problema del semiespacio ante una carga dinámi-ca tipo escalón concentrada, en el presente apartado se muestra la solución presentada

3.5 respuesta dentro del semiespacio 43

por Johnson [19] para el caso más general en el que tanto la fuente como el receptorse encuentran dentro del semiespacio. Es evidente que esta solución es más complejapuesto que requiere realizar la transformada inversa de Laplace L−1(ξx ,ξy,s) a todos lostérminos de la Ecuación (3.25). Sin embargo, el modo de proceder es totalmente análo-go respecto al caso particular, donde el procedimiento de Cagniard-de Hoop descritopreviamente es completamente válido. Habida cuenta de ello, aquí solo se expone elresultado final.

G(x, y, z, t; 0, 0, z′) =[P(x, y, z, t; 0, 0, z′) + S(x, y, z, t; 0, 0, z′) + PP(x, y, z, t; 0, 0, z′)

+SS(x, y, z, t; 0, 0, z′) + PS(x, y, z, t; 0, 0, z′) + SP(x, y, z, t; 0, 0, z′)]B (3.73)

Empezando por las ondas directas P y S, la transformada inversa es sencilla dedesarrollar sin más que aplicar el método. De hecho la integral sobre la variable deintegración p puede evaluarse de forma analítica tal como se expone a continuación.Es evidente que la suma de sendas soluciones corresponde con la solución del espaciocompleto.

Figura 3.13: Representación del semiespacio del problema general relativo a los términos de las ondasdirectas P y S

P(x, y, z, t; 0, 0, z′) =1

8πρrD(x, y, z, t, z′)H

(t− r/cp

)(3.74)

S(x, y, z, t; 0, 0, z′) =1

8πρrE(x, y, z, t, z′)H (t− r/cs) (3.75)

siendo:


D11 = (3(t/r)2 − c−2p ) sin2 θ cos2 ϕ − ((t/r)2 − c−2

p )

D12 = (3(t/r)2 − c−2p ) sin2 θ sin ϕ cos2 ϕ

D13 = −(3(t/r)2 − c−2p ) sin θ cos θ cos ϕ

D21 = (3(t/r)2 − c−2p ) sin2 θ sin ϕ cos2 ϕ

D22 = (3(t/r)2 − c−2p ) sin2 θ sin2 ϕ − ((t/r)2 − c−2

p )

D23 = −(3(t/r)2 − c−2p ) sin θ cos θ sin ϕ

D31 = −(3(t/r)2 − c−2p ) sin θ cos θ cos ϕ

D32 = −(3(t/r)2 − c−2p ) sin θ cos θ sin ϕ

D33 = (3(t/r)2 − c−2p ) cos2 θ − ((t/r)2 − c−2

p )

(3.76)

E11 = −(3(t/r)2 − c−2s ) sin2 θ cos2 ϕ + (t/r)2 + c−2

s

E12 = −(3(t/r)2 − c−2s ) sin2 θ sin ϕ cos2 ϕ

E13 = (3(t/r)2 − c−2s ) sin θ cos θ cos ϕ

E21 = −(3(t/r)2 − c−2s ) sin2 θ sin ϕ cos2 ϕ

E22 = −(3(t/r)2 − c−2s ) sin2 θ sin2 ϕ + (t/r)2 − c−2

s

E23 = (3(t/r)2 − c−2s ) sin θ cos θ sin ϕ

E31 = (3(t/r)2 − c−2s ) sin θ cos θ cos ϕ

E32 = (3(t/r)2 − c−2s ) sin θ cos θ sin ϕ

E33 = (3(t/r)2 − c−2s ) cos2 θ − 2((t/r)2 − c−2

s )

(3.77)

En lo que se refiere a las ondas reflectadas PP y SS la transformada inversa se obtieneen forma integral. Puede observarse que sendas soluciones son respectivamente análo-gas a las obtenidas para los términos M y N del caso particular del receptor sobre lasuperficie libre. Incluso todos los comentarios realizados en dicho apartado sobre losproblemas numéricos y formas de paliarlos son válidos aquí.


Figura 3.14: Representación del semiespacio del problema general relativo a los términos de las ondasreflejadas PP y SS

PP(x, y, z, t; 0, 0, z′) =

=1

2π2ρr′

∫√

(t/r′)2−c−2p

0ℜ

1

√(tr′)2 − c−2

p − p2

I(p, q, z, t, z′)

σ

H(

t− r′c−1p

)

dp

(3.78)

SS(x, y, z, t; 0, 0, z′) =

=1

2π2ρr′

∫ p′2

0ℜ

1

√(tr′)2 − c−2

s − p2

J(p, q, z, t, z′)

σ

H(t− t′2

)dp

(3.79)

Siendo:


¯I11 = −(γ2 − 4ηpηs(q2 − p2))((q2 + p2) cos2 ϕ − p2)

¯I12 = −(γ2 − 4ηpηs(q2 − p2))(q2 + p2) sin ϕ cos ϕ

¯I13 = −(γ2 − 4ηpηs(q2 − p2))qηp cos ϕ

¯I21 = −(γ2 − 4ηpηs(q2 − p2))(q2 + p2) sin ϕ cos ϕ

¯I22 = −(γ2 − 4ηpηs(q2 − p2))((q2 + p2) sin2 ϕ − p2)

¯I23 = −(γ2 − 4ηpηs(q2 − p2))qηp sin ϕ

¯I31 = −(γ2 − 4ηpηs(q2 − p2))qηp cos ϕ

¯I32 = −(γ2 − 4ηpηs(q2 − p2))qηp sin ϕ

¯I33 = (γ2 − 4ηpηs(q2 − p2))η2p

(3.80)

¯J11 = −(γ2 − 4ηpηs(q2 − p2))((q2 + p2) cos2 ϕ − p2)

+c−2s (γ2 − 4ηpηs(q2 + p2)(2 cos2 ϕ − 1))

¯J12 = −(q2 + p2)(γ2 − 4ηpηs(q2 − p2) + 8c−2s ηpηs) sin ϕ cos ϕ

¯J13 = −qηs(γ2 − 4ηpηs(q2 − p2)) cos ϕ

¯J21 = −(q2 + p2)(γ2 − 4ηpηs(q2 − p2) + 8c−2s ηpηs) sin ϕ cos ϕ

¯J22 = −(γ2 − 4ηpηs(q2 − p2))((q2 + p2) sin2 ϕ − p2)

+c−2s (γ2 − 4ηpηs(q2 + p2)(2 sin2 ϕ − 1))

¯J23 = −qηs(γ2 − 4ηpηs(q2 − p2)) sin ϕ

¯J31 = qηs(γ2 − 4ηpηs(q2 − p2)) cos ϕ

¯J32 = qηs(γ2 − 4ηpηs(q2 − p2)) sin ϕ

¯J33 = −(q2 − p2)(γ2 − 4ηpηs(q2 − p2))

(3.81)

p′2 =

√

(t/r′)2 − c−2s sin θ′ ≤ cs/cp

√(

t/r′−√

c−2s −c−2

p cos θ′

sin θ′

)2

− c−2p sin θ′ > cs/cp

(3.82)

t′2 =

r′/cs sin θ′ ≤ cs/cpr′/cp sin θ′ + r

√

c−2s − c−2

p cos θ′ sin θ′ > cs/cp(3.83)

q =

− tr′ sin θ + ir′ cos θ′

√(tr′)2 − c−2

p − p2 para el término PP

− tr′ sin θ + ir′ cos θ′

√(tr′)2 − c−2

s − p2 para el término SS(3.84)

γ = definido en el apartado anterior (3.85)

σ = definido en el apartado anterior (3.86)


Donde ahora el radio vector se define como :

r′ =√

R+ (z+ z′) (3.87)

θ′ = arctang

(R

z+ z′

)

(3.88)

Los últimos dos términos PS y SP representan una onda que al reflejarse sobrela superficie cambia su naturaleza (S ↔ P). De nuevo, aplicando la transformadainversa mediante el método general puede obtenerse las soluciones en forma integral.No obstante, estos términos exhiben una dificultad añadida debido a que el punto dereflexión depende del tiempo de evaluación. Caginard [8] y Triantafyllidis [48] dancuenta de este problema y exponen una técnica para manejar esta dificultad. En lasiguiente sección se muestran unas pautas sencillas para librar esta dificultad y podercomputar numéricamente las integrales.

Figura 3.15: Representación del semiespacio del problema general relativo a los términos de las ondasreflejadas PS y SP

PS(x, y, z, t; 0, 0, z′) =

=1

2π2ρ

∫ p3

0ℜ[

i

R+ q(z′/ηp + z/ηs)

]K(p, q, z, t, z′)

σ

H (t− t3) dp (3.89)

SP(x, y, z, t; 0, 0, z′) =

=1

2π2ρ

∫ p3

0ℜ[

i

R+ q(z/ηp + z′/ηs)

]L(p, q, z, t, z′)

σ

H (t− t3) dp (3.90)


siendo:

K11 = 4γηs((q2 + p2) cos2 ϕ − p2)

K12 = 4γηs(q2 + p2) sin ϕ cos ϕ

K13 = 4qγηsηp cos ϕ

K21 = 4γηs(q2 + p2) sin ϕ cos ϕ

K22 = 4γηs((q2 + p2) sin2 ϕ − p2)

K23 = 4qγηsηp sin ϕ

K31 = 4qγ(q2 − p2) cos ϕ

K32 = 4qγ(q2 − p2) cos ϕ

K33 = 4qγηp(q2 − p2)

(3.91)

¯L11 = 4γηs((q2 + p2) cos2 ϕ − p2)

¯L12 = 4γηs(q2 + p2) sin ϕ cos ϕ

¯L13 = −4qγ(q2 − p2) cos ϕ

¯L21 = 4γηs(q2 + p2) sin ϕ cos ϕ

¯L22 = 4γηs((q2 + p2) sin2 ϕ − p2)

¯L23 = −4qγ(q2 − p2) cos ϕ

¯L31 = −4qγηsηp cos ϕ ¯L32 = −4qγηsηp sin ϕ ¯L33 = 4qγηp(q2 − p2)

(3.92)

Las variables R, ϕ, ηp, ηs,γ y σ son las mismas que las definidas para los términosprecedentes. En cambio, la variable q se define de forma diferente tal que ahora depen-de del punto de reflexión que varia a su vez con el tiempo.

q = −α + iz′

Rp

√

c−2p + p2

1+ (z′/Rp)2(3.93)

Se han introducido las variables Rp y Rs que representan la proyección horizontal dela trayectoria recorrida por las ondas P y S respectivamente. Evidentemente su sumadebe ser igual al valor de R y conocido cualquiera de ellos, para un instante de tiempo,se conoce el punto de reflexión de la onda. Matemáticamente estos términos deben decumplir la relación siguiente:

Rp + Rs = R

(R+ z′2/Rp + z2/Rs)α = t

(3.94)

donde α es un parámetro definido por:


α2 =1

(1+ (z′/Rp)2

)(1+ (z/Rs)2)

[

cs−2(z/Rs)2(1+ (z′/Rp)2

)− c−2

p (z′/Rp)2(1+ (z/Rs)2

)

(z/Rs)2 − (z′/Rp)2+ p2

]

(3.95)

El valor del limite de integración p3 se calcula resolviendo en primer lugar Rp, Rs

y α en función del tiempo para p = 0 en el sistema de ecuaciones (3.94), y luegosustituyendo en la ecuación siguiente:

α2 =c−2p + p23

1+ (z′/Rp)2(3.96)

Mientras que el valor del tiempo de llega de las ondas t3 se deriva de la mismaforma pero estableciendo p3 = 0 en la Ecuación (3.96).

En la siguiente sección se detalla el procedimiento paso a paso para obtener todos loslímites, tanto el de integración, el tiempo de llegada de las ondas y el lugar geométricode los punto de reflexión.

Por último, observar que los términos PS y PS son similares salvo que la profundi-dad de la fuente y el receptor están intercambiados. De hecho físicamente uno es elinverso del otro, por lo que el procedimiento para calcularlos es el mismo sin más queintercambiar la profundidad de la fuente por la del receptor.

3.5.1 Resultados numéricos

A la vista de la solución presentada, el problema general requiere evaluar numérica-mente seis términos relativos a las ondas P,S,PP,SS,SP y PS. Los dos primeros términos(ondas P y S) están deducidos en forma cerrada, y carecen de problemas numéricos.Los siguientes dos términos relativos a las ondas reflejadas PP y SS, desde el punto devista matemático, son semejantes a los términos M y N respectivamente del problemaparticular superficial, sin más que tener en cuenta que ahora el ángulo de emisión sedenota por θ′ y el radio r′. De hecho se puede verificar que comparten los mismo pro-blemas numéricos, cuya forma de resolverlos ya se detalló en el apartado 3.4.1. Aquíse limita a resolver los nuevos problemas que surgen de los términos PS y SP. Fíjeseque estos términos están exentos de singularidad en el integrando, salvo cuando tantofuente como receptor tienden a la superficie y excitan los polos de Rayleigh (la formade proceder fue ya presentada). El problema aparece debido al hecho de que el puntode reflexión no es fijo y depende del tiempo de evaluación, es decir, existe un lugargeométrico de posibles puntos de reflexión tal como se muestra en la Figura 3.15. Estosignifica que el límite de integración (p3) y la variable principal q que dependen del


punto de reflexión, también dependa del tiempo de evaluación de forma indirecta. Ensi mismo esto no supone un problema numérico, pero tal como se vislumbra de la for-mulación presentada y como se verá a continuación, averiguar el punto de reflexiónconlleva resolver un sistema de ecuaciones paramétricas no lineal, cuyos parámetrosson el tiempo de evaluación y la variable de integración p, que requiere de ciertastécnicas numéricas. Cabe mencionar de nuevo que ambos términos PS y SP son unoel inverso del otro, con lo cual comparten los mismos procedimientos. De aquí enadelante se referirá al término PS.Aquí conviene aclarar cierta nomenclatura que se va a utilizar (Figura 3.15). El án-

gulo de emisión, conocido como el ángulo formado por la dirección del frente de ondacon la vertical, que en este caso es variable con el tiempo, se denota con i1, mientrasque el ángulo que incide en el receptor tras la reflexión es i2. Por otro lado la distan-cia horizontal entre fuente y receptor R se divide ahora en la recorrida por las ondasprimarias Rp y la recorrida por las secundaria Rs. para mejor aclaración.

Aclarado esto, parece lógico calcular en primera instancia el lugar geométrico delos puntos posibles de reflexión definidos por los valores límites de Rp o Rs. El pun-to de reflexión en función del tiempo se localizará dentro de esos límites. Sin másque resolver un problema de minimización o maximización de tiempos se obtiene lasrelaciones:

min(t) = min[rp/cp + rs/cs

]

⇒ sin i1 =cp

cssin i2 (3.97)

que determina el limite inferior de Rs o lo que es lo mismo el límite superior de Rp.Realizando la misma operación pero maximizando el tiempo se obtiene evidentementeel límite superior de Rs, y viene definido por la relación:

i1 = i2 (3.98)

Ahora es lógico calcular los limites de integración p3(t) en función del tiempo, y eltiempo de llegada de las ondas t3. Evidentemente él tiempo de llega de las ondas debecoincidir con el menor tiempo hallado en la minimización anterior. Alternativamentepuede ser obtenido resolviendo la ecuación (3.96) con p3 = 0 y p = 0, y determinandolos valores de Rs, Rp, α. Sustituyendo en (3.94) se obtiene el valor de t3. Notar quela ecuación que se resuelve es muy compleja y no lineal, lo que requiere técnicasnuméricas. En este caso, estos cálculos pueden evitarse al calcularse los límites deRp y Rs, sin embargo no puede hacerse lo mismo con el límite de integración p3que depende del tiempo de evaluación. Nuevamente hay que resolver el sistema deEcuaciones (3.94) con p = 0 y obtener Rs, Rp, α en función del tiempo. Sustituyendoen la Ecuación (3.96) se obtiene el valor del límite de integración. Notar que el sistema


puede ser reducido a una sola ecuación sin más que utilizar la relación R = Rp + Rs

y dejando una sola incógnita con el tiempo como parámetro (por ejemplo Rs = f (t)),pero es fuertemente no lineal con múltiples raíces. En este sentido aquí se planteandos procedimientos para resolver este problema:

Desde el punto de vista puramente numérico puede resolverse utilizando el Mé-todo de Newton Raphson. Es un método iterativo que linealiza la ecuación encada punto de colocación relativo a cada iteración hasta llegar un valor de conver-gencia deseado. Este método es sencillo y robusto, pero inestable cuando existenpuntos límites. Aquí se propone el Método de Arc-Length capaz de pasar pun-tos límites. La diferencia radica que en cada iteración el avance hacia la raíz estadetermino por un arco. Para más detalle véase Crisfield1 [10].

Desde el punto de vista analítico las raíces de la ecuación pueden reducirse alas raíces de una ecuación polinómica. Evidentemente éstas raíces pueden serhalladas más fácilmente, e incluso puede plantearse como un simple problemade autovalores. La raíz buscada es aquella que se encuentra dentro de los límitesestablecidos anteriormente de Rs y Rp. La ecuación en forma paramétrica parael término PS es (para el término SP basta con cambiar z ↔ z′):

GR6s + FR5

s + ER4s + DR3

s + CR2s + BRs + A = 0

(3.99)

siendo:


A = c2pc2sR4t

2z4 + (−1)c2pR4z′6 + c2pc

2sR

2t2z′4z′2 + (−1)c2pR2z6z′2 + c2sR

2z6z′2

B = ((−2)c2pR5z4 + (−4)c2pc

2sR

3t2z4 + 4c2pR3z6 + (−4)c2pR

3z4z′2 + 2c2sR3

z4z′2 + (−2)c2pc2sRt

2z4z′2 + 2c2pRz6z′2 + (−2)c2sRz

6z′2 + (−2)

c2pRz4z′4 + 2c2sRz

4z′4)

C = ((−1)c2pR6z2 + c2pc

2sR

4t2z2 + 8c2pR4z4 + 6c2pc

2sR

2t2z4 + (−6)c2pR2z6 + (

− 3)c2pR4z2z′2 + c2sR

4z2z′2 + 10c2pR2z4z′2 + (−3)c2sR

2z4

z′2 + c2pc2s t

2z4z′2 + (−1)c2pz6z′2 + c2s z

6z′2 + (−3)c2pR2

z2z′4 + 2c2sR2z2z′4 + (−1)c2pc

2s t

2z2z′4 + 2c2pz4z′4 + (−2)

c2sz4z′4 + (−1)c2pz

2z′6 + c2s z2z′6)

D = (4c2pR5z2 + (−4)c2pc

2sR

3t2z2 + (−12)c2pR3z4 + (−4)c2pc

2sRt

2z4 + 4c2pR

z6 + 8c2pR3z2z′2 + (−8)c2pRz

4z′2 + 4c2pRz2z′4)

E = ((−6)c2pR4z2 + 6c2pc

2sR

2t2z2 + 8c2pR2z4 + c2pc

2s t

2z4 + (−1)

c2pz6 + c2sR

4z′2 + (−1)c2pc2sR

2t2z′2 + (−7)c2pR2z2z′2 + (−3)c2s

R2z2z′2 + 2c2pz4z′2 + c2s z

4z′2 + 2c2sR2z′4 + (−1)c2pc

2s t

2

z′4 + (−1)c2pz2z′4 + (−2)c2s z

2z′4 + c2sz′6)

F = (4c2pR3z2 + (−4)c2pc

2sRt

2z2 + (−2)c2pRz4 + (−2)c2sR

3z′2 + 2c2pc2sRt

2z′2 + 2c2pRz2z′2 + 2c2sRz

2

z′2 + (−2)c2sRz′4)

G = ((−1)c2pR2z2 + c2pc

2s t

2z2 + c2sR2z′2 + (−1)c2pc

2s t

2z′2) (3.100)

Por último notar que la variable q(p, t) depende del tiempo y de la variable de inte-gración p indirectamente a través del punto de reflexión. De nuevo hay que resolverel sistema de ecuaciones (3.94) y obtener Rs(t, p), Rp(t, p), α(t, p). La forma de proce-der es exactamente igual que para obtener el límite de integración p3 pero ahora hayque añadirle el parámetro p a parte del parámetro temporal t a la ecuación a resol-ver. Igualmente se puede resolver de forma numérica o analítica, pero mediante éstaúltima se llegan que se recogen en el anexo A. En la Figura 3.16 se muestra una compa-rativa típica de ambos procedimientos para averiguar el valor de q. Puede observarseque se alcanza la misma solución, pero el procedimiento analítico tarda tan solo el1.4% del tiempo que necesitado el numérico. Sin embargo la solución analítica es muyoscilatoria e inestable en ciertos casos.En las Figuras 3.17 y 3.18 se muestra la evolución temporal de los desplazamiento

de un semiespacio con las mismas propiedades dinámicas que las tomadas en el apar-tado 3.4.1. Corresponde a las dos mismos casos que los reproducidos por Johnson [19]


0 1 2 3 4 5 6

x 10−4

4.31

4.312

4.314

4.316

4.318

4.32

4.322

4.324

4.326

4.328x 10

−4

p [−]

q(p)

[−]

M. numericoM. analitico

Figura 3.16: Evaluación de la variable q(p, t) mediante procedimiento numérico y analítico.

y Triantafyllidis [48], donde tanto la fuente como el receptor se localizan dentro del se-miespacio. A parte de las observaciones realizadas anteriormente, en la Figura 3.17 a)se puede ver la complejidad de la solución general dentro del semiespacio, en dondese señalan el tiempo de llegada de todas las ondas. Con esta solución se puede repre-sentar todos los posibles casos y se puede ver cómo los términos se cancelan cuandotiende la fuente y/o el receptor a la superficie. Sin embargo, evidentemente ésta es mu-cho más costosa y compleja numéricamente que la solución particular anteriormentepresentada.


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

SSS

SPPPP

PS

(a) GHxx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]D

espl

azam

ient

o [−

]

(b) GHxz

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) GHyy

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) GHzx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(e) GHzz

Figura 3.17: Evolución temporal de las componentes de la función Green adimensional, GH(2× 103, 0, 1×103, t; 0, 0, 2× 103), de un semiespacio homogéneo.


0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

SSSP

PPSP

PS

(a) GHxx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) GHxz

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) GHyy

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) GHzx

0 0.5 1 1.5 2

−0.5

0

0.5

Tiempo [−]

Des

plaz

amie

nto

[−]

(e) GHzz

Figura 3.18: Evolución temporal de las componentes de la función Green adimensional, GH(2× 103, 0, 4×103 , t; 0, 0, 2× 103), de un semiespacio homogéneo.

4SOLUCIÓN FUNDAMENTAL DEL SEMIESPACIOESTRAT IF ICADO EN EL DOMIN IO DEL T IEMPO

4.1 introducción

En el capítulo anterior se ha presentado la respuesta de un semiespacio homogéneoante una carga impulsiva. Se ha demostrado que aún siendo un problema relativamen-te sencillo, la solución es realmente compleja, presentándose en forma integral con nu-merosas dificultades. Para el mismo problema dinámico pero para un semiespacio ho-rizontalmente estratificado no existe solución analítica en el dominio espacio-tiempo.Pueden encontrarse numerosos artículos que han tratado de resolver este problema,empezando por el trabajo desarrollado en 1950 por Thomson y Haskell [47] quienesintrodujeron el método de la matriz de propagación o de transferencia. Este métodopuede considerarse como el precursor, y se basa en realizar una discretización del me-dio según los diferentes estratos que lo componen, de forma que se monta una matrizde propagación que compatibiliza las velocidades y tensiones en las interfases conti-guas. Posteriormente han surgido nuevos métodos, entre ellos destacar el método dela transformación integral que consiste en elevar el problema original en el espacio-tiempo al dominio de número de onda-frecuencia, donde las ecuaciones simplificansustancialmente. En esta línea citar los trabajos desarrollados por Buouchon [7], Lu-co and Aspel [3, 4], Xu and Mal [49] y Hisada [17, 18]. Basado en este método peroformulado de forma sistemática debe mencionarse el Método de la Matriz de Rigidezdesarrollado por Kausel y Roëset [26] en 1981. En el dominio transformado númerode onda-frecuencia existen soluciones a las ecuaciones de Navier tanto para el semies-pacio homogéneo como para una placa estratificada. El método de Kausel y Roëset lasutiliza para plantear matrices elementales que relacionan los desplazamientos en lasinterfases con las fuerzas aplicadas sobre ellas. La matriz de rigidez global del semies-pacio estratificado horizontalmente se obtiene tras ensamblar las matrices elementalessin más que aplicar compatibilidad de desplazamientos entre estratos. Finalmente sellega a una estructura P=KU donde P y U son las tracciones y desplazamientos enlas interfases respectivamente, y cuya solución es exacta. Este método requiere de unatransformación inversa del dominio número de onda-frecuencia a espacio-tiempo, quepresenta numerosas dificultades. Alternativamente, Kausel (1994) [22] propuso una va-riante discreta de este método llamada "Thin-Layer Method"(TLM). Parte de la mismaidea pero el medio se discretiza en delgada láminas en el sentido de la estratificación,y en el que ahora se propone funciones de forma para aproximar los desplazamientos.La ventaja principal es que se llega a matrices más sencillas, lo que permite plantear

57

58 solución fundamental del semiespacio estratificado en el dominio del tiempo

la solución mediante el método de superposición modal tras resolver un problemaestándar de autovalores y autovectores. Es decir, directamente se obtiene en el númerode onda-tiempo, y donde la transformada inversa al espacio carece de singularida-des. La desventaja de este método es que a priori no es aplicable a un semiespacioestratificado sino sólo a una placa horizontalmente estratificado con espesor finito. Re-cientemente Park y Kausel (2006) [24, 25] han resuelto este problema al proponer unmétodo híbrido que combina la solución de la placa horizontalmente estratificada conel semiespacio homogéneo en el dominio del número de onda-tiempo.En este capítulo se expone la formulación del método híbrido propuesto por Park y

Kausel [24, 25] relativa a la función de Green del semiespacio estratificado tridimensio-nal ante una carga impulsiva. Este método se apoya en dos pilares maestros que son:la función de Green del semiespacio homogéneo en el número de onda-tiempo (deforma abreviada k-t), que se conoce analíticamente [37]; y la función de Green de laplaca horizontalmente estratificada en el número de onda-tiempo obtenida medianteel TLM [22]. Ambas soluciones se desgranan en los siguientes dos subcapítulos. Porotro lado, tal como se comentó en el capítulo 2, la solución tridimensional en coorde-nadas cilíndricas puede plantearse como una combinación adecuada de las solucionesal problema de deformación plana. De hecho, en éste último se da la circunstanciaque los desplazamientos en el plano se desacoplan de los desplazamientos fuera de él,dividiéndose este problema a su vez en dos relativos a las ondas P-SV y SH respectiva-mente. Una vez conocidas las dos soluciones maestras, en el subcapítulo 4.4 se exponeel acoplamiento entre ambas, dando lugar a las funciones de Green bidimensionales(P-SV y SH) del semiespacio estratificado en el número de onda-tiempo (k-t) [24, 25].La función de Green tridimensional en el dominio espacio-tiempo se calcula finalmen-te mediante combinación de las funciones de Green bidimensionales y realizando unatransformada inversa de Hankel [43]. Como se deja entrever, este método no es incon-dicionalmente estable: en parte se debe al acoplamiento de las dos soluciones maestrasque se computan de formas diferentes; y en parte debido al TLM que requiere de unacierta discretización vertical espacial de los estratos. Por ello, él último subcapítulo sededica a estudiar la precisión y estabilidad del método.Aquí se cree conveniente mostrar un esquema de la resolución que se sigue. La idea

es mostrar al lector una visión global del método y elaborar una guía del procedimien-to a desarrollar. Básicamente el método consiste en los siguientes pasos:

1. Calcular las funciones de Green en el dominio k-t del semiespacio homogéneo2D ante una carga concentrada impulsiva para el problema en el plano (P-SV)y fuera del él (SH) [37]. La solución a ésta última se conoce de forma cerrada,mientras que para las ondas P-SV debe plantearse una transformada inversa deFourier sobre la frecuencia. Esta integral presenta singularidades debido a lospolos de Rayleigh, que se evitan tomando un camino complejo de integración, yademás son fuertemente oscilatorias. Pese a la complejidad y alto coste compu-

4.2 solución del semiespacio homogéneo 59

tacional que requieren pueden pre-computarse y almacenarse en tablas para suposterior uso.

2. Calcular las funciones de Green P-SV y SH en el dominio k-t de una placa ho-rizontalmente estratificado ante una carga impulsiva [35]. Ambas soluciones secalculan con el TLM en el dominio temporal. Se construye la matriz de rigi-dez de la placa en dirección normal a la superficie y en el dominio número deonda-frecuencia. Se calcula las frecuencias y modos de vibración y se plantea lasolución en el dominio temporal mediante el método de la superposición modal.

3. Acoplamiento de la función de Green del semiespacio con la de la placa estra-tificada relativo al problema P-SV y SH [24, 25]. El acoplamiento se realiza re-solviendo un sistema de ecuaciones en forma de convolución temporal. De estaforma se alcanzan las funciones de Green P-SV y SH de un semiespacio hori-zontalmente estratificado. La región estratificada debe ser finita por requisito delTLM, mientras que el resto se considera homogéneo.

4. Cálculo de la función de Green ante una carga puntual en el dominio númerode onda radial y tiempo (Kr-t) por combinación de las funciones de Green bidi-mensionales P-SV y SH [43] . Por último se realiza una Transformada inversa deHankel.

4.2 solución del semiespacio homogéneo

En este apartado se muestran las funciones de Green del semiespacio homogéneopara los problemas P-SV y SH en el dominio número de onda horizontal-tiempo (k-t),y que fue obtenida por Park y Kausel (2004) [37]. Estas funciones tienen gran inte-rés puesto que son aplicables para cualquier tipo de carga con distribución aleatoriahorizontal, y por supuesto es válida para ser aplicada al método híbrido.

En términos generales para este apartado, se considera un semiespacio bidimensio-nal (plano x-z) definido en la región negativa del eje z y cuyo vector director ~ezesnormal a la superficie en dirección vertical. Notar que se ha tomado un vector directororientado hacia fuera del semiespacio, en contra de la notación tomada en el aparta-do (3). Esta convención se toma en aras a facilitar la combinación de la solución delsemiespacio con la placa realizada posteriormente, en donde se toma como referenciala interfase de unión. Véase Figura 4.1

4.2.1 Problema bidimensional en el plano (SH)

Como se expuso en el capítulo 2, la ecuación que gobierna la respuesta de un semi-espacio en dirección fuera del plano uy ante una carga SH temporal y horizontalmentearmónica aplicada en la superficie, es la siguiente:


Figura 4.1: Geometría del semiespacio bidimensional: Superficie libre coincide con el plano z = 0 tal queel eje z negativo define la región del semiespacio.

∂2uy

∂x2+

∂2uy

∂z2=

1c2s

∂2uy

∂t2(4.1)

donde se ha supuesto que la fuerza se introduce como condición de contorno y nocomo fuerza de volumen.La solución bien conocida [23, 37], es una función horizontal y temporalmente ar-

mónica que decae exponencialmente con la profundidad, y tal como se demostró enel mencionado capítulo, puede expresarse genéricamente por:

uy = ( ISHe−ikzs + RSHe

+ikzs)ei(ωt−kxx) (4.2)

Donde los coeficientes ISH, RSH son constantes de integración aún desconocidas, ykzs es el número de onda vertical que se relaciona con el número de onda horizontalsegún la expresión siguiente obtenida tras sustituir en la ecuación (4.1):

ksz =

√

ω2

c2s− k2x =

√

k2s − k2x = kx

√√√√−

(

1−(kskx

)2)

= ikxs (4.3)

Teniendo en cuenta que el semiespacio está definido en la región negativa del eje z,y se propaga en la misma dirección, la respuesta se reduce a:

uy = ( ISHe−ikzs)ei(ωt−kxx) = ( ISHe

zkxs)ei(ωt−kxx) (4.4)

Que evidentemente satisface la condición de radiación para z → ∞. Por otro lado, latensión tangencial en dirección ~ey para un plano horizontal de normal ~ez viene dadapor:


τyz = µ∂uy

∂z= µkxsuy (4.5)

La constante de integración puede deducirse sin más que aplicar la condición decontorno en tensión sobre la superficie del semiespacio:

Py(kx,ω) = Pei(ωt−kxx) = µksISHezkxs

=⇒ ISH =P

µkxs(4.6)

Sustituyendo en la Ecuación 4.4, el campo de desplazamientos resultante es:

uy == (P

µkxsezkxs)ei(ωt−kxx) = PGδ

yyei(ωt−kxx) (4.7)

Y por tanto la función de Green en el dominio k− ω para una carga de magnitudunidad puede expresarse como:

Gyy(kx,ω) =ezkxs

µkxs=

ez√

k2x−k2s

µ√

k2x − k2s(4.8)

Realizando ahora una transformada inversa de Fourier sobre la frecuencia, se obtie-ne finalmente la función de Green en el dominio k− t buscada:

Gδyy(kx, t) =

∫ ∞

∞

ezkxs

µkxs=

1csρ

J0

(

kxcs

√

t2 − t2s

)

H(t− ts) (4.9)

Donde J0 es la función de Bessel de orden 0 y ts = |z|/cs es el tiempo de llega de lasondas S. El símbolo ˇ sobre la variable denota la trasformada de Fourier al dominionúmero de onda-tiempo.

Notar que la función de Green hallada es válida para un semiespacio elástico sinatenuación, no obstante es posible modificar la solución para introducirle amortigua-miento [37]. Existen diferente tipos de modelos, el más sencillo de introducir es elamortiguamiento viscoso proporcional a la masa y dependiente del número de onda.Éste, bajo cierta condiciones o más bien elecciones, puede reproducir además otrosmodelos, por ejemplo, amortiguamiento histéretico o amortiguamiento material. Enprincipio, físicamente no es posible incorporar amortiguamiento independiente de lafrecuencia sobre la formulación directamente en el tiempo, puesto que viola el prin-cipio de causalidad, sin embargo puede simularse con una excelente aproximación, ysiendo despreciable la inexactitud física de la solución. Para ello se parte de la misma


ecuación del movimiento 4.1 pero se le añade un término inelástico que representa elamortiguamiento viscoso proporcional a la masa. La ecuación queda entonces de laforma:

∂2uy

∂x2+

∂2uy

∂z2=

1c2s

∂2uy

∂t2+ D

∂uy

∂t(4.10)

Donde el coeficiente de amortiguamiento puede definirse por D = 2ρη, siendo η uncoeficiente de dimensión frecuencial.La solución a la ecuación es análoga al caso elástico anterior, pero ahora el número

de onda vertical se ve afectado por el coeficiente de amortiguamiento según la expre-sión:

k′sz =

√√√√√√

(ω − iη

cs

)2

︸︷︷︸

k′2s

−(

k2x −η2

c2s

)

︸︷︷︸

k′2x

(4.11)

Que reordenando y renombrado las variables puede expresarse también de formaanáloga al caso elástico:

k′sz =

√(

ω′

cs

)2

− k′2x (4.12)

Notar que ahora la frecuencia se convierte compleja, en el que se añade una constan-te imaginaria iη, y el número de onda se desplaza por la constante η2/c2s . Por últimorealizando una transformada inversa sobre la frecuencia en el plano complejo, dondese ha realizado el cambio de variable ω → ω′, se obtiene la función de Green en eldominio k − t. Ésta es análoga al caso elástico, pero ahora aparece un nuevo factorexponencial e−ηt y el número de onda se encuentra desplazado:

Gδyy(kx, t) = e−ηt 1

csρJ0

(

kxcs

√

t2 − t2s

)

H(t− ts) (4.13)

Si se escoge η = ξkxcs, siendo ξ la fracción del amortiguamiento crítico, y sustitu-yendo k′x =

√

k2x − (η/cs)2, se llega finalmente a la expresión exacta de la función deGreen de un semiespacio homogéneo inelástico con amortiguamiento viscoso propor-cional a la masa.

Gδyy(kx, t) = e−ξkxcst

1csρ

J0

(

kxcs

√

t2 − t2s

)

H(t− ts) (4.14)


Efectivamente observar que ξ no depende de la frecuencia, mientras que puede de-pender del número de onda. Si se toma adicionalmente constante e independiente deéste último, efectivamente se puede reproducir el amortiguamiento material indepen-diente de la frecuencia.

En la Figura 4.2 se muestra la evolución temporal de la función de Green SH anteuna carga impulsiva tipo senoidal. Se ha considerado un semiespacio cuyas propida-des dinámicas son: cp = 467.70 m/s, cs = 250 m/s, ρ = 1850 kg/m3 y sin amorti-guamiento. La fuente y el receptor son superficiales,y distan 10 m. Se ha utilizado laformulación desarrollada en este apartado y se ha realizado la transformada inversade Fourier sobre el número de onda numéricamente de forma discreta. Además seha superpuesto la función de Green SH en el dominio natural espacio-tiempo exactapropuesta por Eringen y Suhubi [16] para un impulsivo delta de Dirac. Notar que larespuesta se presenta en forma adimensional Gδ

yy = Gyyµπr/Z, siendo Z la magnitudde la fuerza, frente al tiempo adimensional τ = tcs/r. Puede apreciarse que ambas so-luciones prácticamente coinciden, y debido a la adimensionalización las ondas S, queson las únicas que pueden registrarse, llega para el tiempo correcto τ = 1

0 0.5 1 1.5 2−1000

−500

0

500

1000

Tiempo [−]

Des

plaz

amie

nto

[−]

(a) Gδyy

Figura 4.2: Evolución temporal de las componentes de la función Green adimensional, Gδ(10, 0, 0, t; 0, 0, 0)en el plano (P-SV) de un semiespacio homogéneo. Linea continua negra corresponde con la función deGreen propuesta por Park y Kausel, y los círculos grises corresponde con la exacta propuesta por Eringeny Suhubi

4.2.2 Problema bidimensional anti-plano (P-SV)

Aunque la forma de proceder para obtener la función de Green del semiespacio P-SV en el dominio k-t es similar al problema fuera del plano, la complejidad no es nipor atisbo comparable. Igualmente se plantea el problema en el domino transformadonúmero de onda-frecuencia, cuya solución tal como Kausel [37] demuestra es sencillay cerrada. La respuesta en el dominio temporal se alcanza nuevamente mediante unaantitransformada de Fourier sobre la frecuencia. Aquí radica precisamente la comple-


jidad de la solución puesto que en el integrando aparecen singularidades asociadasa las ondas de Rayleigh, que obliga a utilizar técnicas matemáticas para evitarlas. Eneste apartado se muestra el desarrollo de la solución propuesta por Park y Kausel[37, 25] , que tal como se comentó en la introducción, se presenta en forma integral ydebe ser evaluada numéricamente. Aún careciendo de singularidades, en esta secciónse demuestra que su evaluación es dificultosa, debido a la naturaleza oscilatoria delnúcleo y a que es impropia, de ahí que propongan algunas técnicas para paliar estosproblemas.

4.2.2.1 Función de Green en el dominio número de onda-frecuencia

Considerando un semiespacio elástico en ausencia de fuerzas de volumen, las ecua-ciones que gobiernan el campo de desplazamientos, tal como se expresó en el capitulo2, son las siguientes:

[

(λ + 2µ) ∂2

∂x2+ µ ∂2

∂z2(λ + µ) ∂2

∂x∂z

(λ + µ) ∂2

∂x∂z µ ∂2

∂x2+ (λ + 2µ) ∂2

∂z2

] [

ux

uy

]

= ρ

[

¨xu¨yu

]

(4.15)

Se supone que la fuerza se aplica como condición de contorno sobre la superficieen cualquiera de las direcciones del plano, y es temporal y horizontalmente armónica.Por lo tanto la solución bien conocida que satisface la condición de radiación es:

[

ux

uy

]

=

[

1 −is

ip 1

] [

ekx pz 1

1 ekxsz

] [

Ix(P−SV)

Iz(P−SV)

]

ei(ωt−kxx) (4.16)

donde se ha tenido en cuenta que el número de onda vertical de las ondas P y SVrespectivamente satisfacen la relación:

ksz =√

k2s − k2x = kx

√√√√−

(

1−(kskx

)2)

= ikxs

kpz =√

k2p − k2x = kx

√√√√−

(

1−(kp

kx

)2)

= ikxp (4.17)

s =

√√√√

(

1−(kskx

)2)

p =

√√√√

(

1−(kp

kx

)2)

(4.18)

Nuevamente el vector IP−SV son una pareja de constantes de integración y derivan deimponer las condiciones de contorno en tensión sobre la superficie, esto es:


[

Px

Pz

]

=

[

µ ∂∂z µ ∂

∂x

λ ∂∂x (λ + 2µ) ∂

∂z

] [

ux

uz

]

z=0

(4.19)

Reagrupando términos y eliminando el término exponencial ei(ωt−kxx) , se deduce lafunción de Green ante una carga impulsiva en el dominio k− ω [37]. Como se introdu-jo, es un tensor bidimensional que se expresa en forma cerrada por las componentes:

Gxx(kx,ω) =1

4kxµ∆(2sekx pz − s(1+ s2)ekxsz)

Gzx(kx,ω) =i

4kxµ∆(2pekx pz − (1+ s2)ekxsz)

Gxz(kx,ω) =i

4kxµ∆((1+ s2)ekx pz − 2psekx sz)

Gzz(kx,ω) =1

4kxµ∆(−p(1+ s2)ekx pz + 2pekx sz) (4.20)

Donde ∆ es la conocida función de Rayleigh tan problemática, y el número de ondavertical relativo a las ondas P y SV deben cumplir las siguientes condiciones parasatisfacer la condición de radiación:

∆ = ps− 1/4(1+ s2)2 (4.21)

ℜ[kx p] ≥ 0

ℑ[kx p] ≥ 0

ℜ[kxs] ≥ 0

ℑ[kxs] ≥ 0

(4.22)

4.2.2.2 Función de Green en el dominio número de onda-tiempo

Formalmente, la función de Green del semiespacio homogéneo P-SV en el dominiok− t se alcanza al realizar la trasformada inversa de Fourier de la función de Green enel dominio k−ω sobre la frecuencia [37, 35]. Es conveniente introducir aquí la notaciónadimensional puesto que en primer lugar simplifica propiamente la notación, y ensegundo lugar, y más importante, hace extensible la función de Green calculada paraun valor determinado de las variables adimensionales a amplio rango de semiespacioscon diferentes propiedades dinámicas. Dicho esto se define:


a = cs/cp =

√

1− 2ν

2− 2νrelación entre las ondas S y P

γ = cR/cs relación de las ondas de Rayleigh y S

τ = kxcst tiempo adimensional

τp = kxa|z| tiempo adimensional de llegada de las ondas P

τs = kx|z| tiempo adimensional de llegada de las ondas S

Ω =ω

kxcsfrecuencia adimensional

s =√

1− Ω2 número de onda vertical de las ondas P adimensional

p =√

1− a2Ω2 número de onda vertical de las ondas S adimensional

∆ = sp− (1− 12

Ω2)2 Función de Rayleigh adimensional

(4.23)

Luego la transformada inversa de Fourier de la función Green sobre el número deonda, en términos adimensionales puede ser escrita de la forma:

Gδij(kx, z, t) =

18πρcs

∫ ∞

∞Kij(s, p,∆)eiΩτdΩ, i, j = x, z (4.24)

donde los núcleos correspondientes son:

Kxx(kx,ω) =1∆(2sekx pz − s(1+ s2)ekxsz)

Kzx(kx,ω) =i

∆(2pekx pz − (1+ s2)ekxsz)

Kxz(kx,ω) =i

µ∆((1+ s2)ekx pz − 2psekx sz)

Kzz(kx,ω) =1

µ∆(−p(1+ s2)ekx pz + 2pekxsz) (4.25)

Los núcleos de integración son singulares debido los polos reales de la función deRayleigh. Para evaluar las integrales es necesario realizar una integración en el planocomplejo utilizando el Teorema de los Residuos. En este sentido Park y Kausel [37]propusieron tomar el camino de integración mostrado en la Figura 4.3 que consiste en:el eje real de la frecuencia adimensional, que es justamente el camino de integraciónque se pretende hallar; regreso bajo el eje real hasta el punto ΩS, semicírculo alrededorde éste último; recorrido sobre el eje real hasta el infinito; curva compleja de cierrealrededor del infinito; y camino simétrico en la parte negativa del eje de abscisas. Notar


que el camino se ha deformado alrededor de los puntos ΩS y ΩP que correspondecon los ceros de los números de ondas verticales (s(Ω) = 0, p(Ω) = 0) y cuyosvalores son ΩS = 1 y ΩP = 1/a. El motivo es evitar los puntos no analíticos de losnúcleos, condición necesaria para poder aplicar el Teorema de los Residuos. Por otrolado observar también que la curva se ha cerrado alrededor del plano complejo cuyaparte imaginaria es positiva para satisfacer la condición de radiación. Por tanto, laintegral deseada a lo largo del eje real se puede calcular como la suma de los residuosen los polos de Rayleigh menos la integral a lo largo del camino restante en el sentidode las agujas del reloj:

Figura 4.3: Contorno complejo de integración de la función de Green. Park y Kausel [37]

Gδij(kx, z, t) = 2πi ∑ Residuo

[Gij(kx, z,Ω),ΩR

]−

0∫

Ωc,∞

Gij(kx, z,Ω)dΩ

Rama integral︷︸︸︷

−∫ 1/a

∞Gij(kx, z,Ω)dΩ −

∫ 1

1/aGij(kx, z,Ω)dΩ −

∫ 1/a

1Gij(kx, z,Ω)dΩ

−∫ ∞

1/aGij(kx, z,Ω)dΩ −

∫ −1/a

−∞Gij(kx, z,Ω)dΩ −

∫ −1

−1/aGij(kx, z,Ω)dΩ

−∫ −1/a

−1Gij(kx, z,Ω)dΩ −

∫ −∞

−1/aGij(kx, z,Ω)dΩ (4.26)

Se puede demostrar que la integral de cierre sobre el plano complejo es cero, eigual con las integrales alrededor de los puntos ΩS y ΩP (Kausel [37]). Por tanto latransformada inversa de Fourier puede reducirse a:

Gδij(kx, z, t) =

1ρcs

[Rij − Bij] (4.27)

En el que Rij contiene los residuos asociados a los polos de Rayleigh, y Bij se refiere ala rama integral. A continuación se expresa el resultado final de ambos términos, sudesarrollo más minucioso puede encontrarse en Park y Kausel [37].


residuo :

Rij(γ, τ, d) =Nij(γ, d)D(γ)

sinγτ i.j = x, z (4.28)

Siendo:

Nxx(γ, d) = sR(e−τspR − qRe

−τssR)

Nzx(γ, d) = i(pRsRe−τspR − qRe

−τssR)

Nxz(γ, d) = i(qRe−τspR − pRsRe

−τssR)

Nzz(γ, d) = pR(−qRe−τspR + e−τssR)

D(γ) = γ

(

1+ a2 − 2a2γ2√

(1− γ2)(1− a2γ2)− (2− γ2)

)

(4.29)

y donde se han definido las variables:

d = |z| profundidad del receptor

sR =√

1− γ2

pR =√

1− a2γ2

número de onda vertical del polo de Rayleigh

qR = 1− 1/2γ2 (4.30)

Fíjese que para el caso particular de z = 0, el numerador Nij se simplifica, por loque el residuo Rij queda solo en función de τ y γ:

Nxx(γ) = 1/2γ2sR

Nzx(γ) = i(pRsR − qR)

Nxz(γ) = −Nzx(γ)

Nzz(γ) = 1/2γ2pR

(4.31)


rama integral:

Bxx =1π

∫ 1/a

1

s

|∆1|2ℜ[∆1(qe

ikx sd − e−kx pd)] sin ΩτdΩ

− 1π

∫ ∞

1/a

s

∆2[cos kx pd− q cos kx sd] sinΩτdΩ

Bxz =i

π

∫ 1/a

1

ps

|∆1|2[qe−kx pd +ℜ(∆1e

ikx sd)] sin ΩτdΩ

+i

π

∫ ∞

1/a

1∆2

[q sin kx pd+ ps sin kx sd] sin ΩτdΩ

Bzx = − i

π

∫ 1/a

1

ps

|∆1|2[ℑ(∆1e

ikx sd)− psqe−kx pd] sin ΩτdΩ

− i

π

∫ ∞

1/a

1∆2

[ps sin kx pd+ q sin kx sd] sin ΩτdΩ

Bzz =1π

∫ 1/a

1

p

|∆1|2ℑ[∆1(e

ikx sd − qe−kx pd)] sin ΩτdΩ

− 1π

∫ ∞

1/a

p

∆2[cos kx sd− q cos kx pd] sinΩτdΩ (4.32)

En el que se han definido las variantes del número de onda y función de Rayleigh:

p =√

a2Ω2 − 1

s =√

Ω2 − 1

∆ = sp− q2

∆1 = −q2 + isp

∆2 = −q2 − is p

q = 1− 1/2Ω (4.33)

Igual, para el caso particular de z = 0 las integrales se simplifican considerablementetal que la rama integral solo depende de τ y γ:

Bxx =12π

∫ 1/a

1

Ω2 sq2

|∆1|2sinΩτdΩ − 1

2π

∫ ∞

1/a

Ω2s

∆2sinΩτdΩ

Bxz =i

2π

∫ 1/a

1

Ω2qps

|∆1|2sinΩτdΩ

Bzx = −Bxz

Bzz =12π

∫ 1/a

1

Ωsp

|∆1|2sinΩτdΩ − 1

2π

∫ ∞

1/a

Ω p

∆2sinΩτdΩ (4.34)


La función de Green mostrada aquí es válida para un semiespacio elástico en ausen-cia de atenuación. Igualmente que para el problema anti-plano SH, es posible introdu-cirle un amortiguamiento a la solución [37]. El procedimiento e interpretación física estotalmente análoga, por lo que en definitiva, en términos matemáticos solo es necesa-rio introducirle un factor exponencial dependiente del la fracción de amortiguamientocrítico que puede depender del número de onda:

Gδij(kx, z, t) =

e−ξτ

ρcs[Rij − Bij] (4.35)

Por otro lado, se debe subrayar dos observaciones sobre la función de Green en eldominio k-t: las componentes Gδ

xx y Gδzz son reales y puramente simétricas respecto va-

lores negativos y positivos del número de onda horizontal, mientras que Gδxz y Gδ

zx sonpuramente imaginarias y antisimétricas, lo cual asegura que la transformada inversade Fourier sobre el número de onda horizontal exista y sea real; según el teorema dereciprocidad la función de Green expuesta es aplicable para el caso inverso en el queel receptor sea superficial y la fuente ubicada dentro del semiespacio

4.2.2.3 Evaluación numérica

En definitiva la función de Green del semiespacio P-SV se compone de dos térmi-nos, uno asociado al residuo en los polos de Rayleigth, y otro que corresponde con larama integral. El primer término evidentemente no muestra dificultad numérica, sinembargo en el término integral se pueden intuir dos problemas: la existencia de unaintegral impropia y la naturaleza oscilatoria del núcleo debido a la presencia de fun-ciones senoidales. De hecho el periodo de oscilación disminuye a medida que aumentael tiempo adimensional (τ) de evaluación. Por lo que ambos problemas suponen unfuerte inconveniente a la hora de computar numéricamente la rama integral.En la Figura 4.4 se muestra el aspecto de los núcleos correspondientes a las inte-

grales Bxx, Bxz, Bzz (Bxz es similar a Bxz) para diferentes valores crecientes del tiempoadimensional y para el caso concreto de z = 0, ν = 0.25. Esta gráfica es totalmenteanáloga para otros valores de z y ν, y por ende es representativa del comportamiento.Los puntos marcados con círculos corresponde a los valores de los núcleos evaluadosen Ωs = 1 y Ωp = 1/a. De la misma se puede extraer dos ideas principales: inde-pendientemente del del tipo de núcleo, exhiben un comportamiento muy diferentespara valores de Ω comprendidos entre [ΩS,ΩP] y [ΩP,∞]; e igualmente, en todos losnúcleos se aprecia una tendencia asintótica a medida Ω tiende a infinito. Ésta últimaobservación puede verse más claramente en la Figura 4.5 donde se representan de nue-vo el comportamiento de los núcleos pero obviando la parte senoidal. En las gráficasse han dibujado también las asíntotas respectivas.La dos observaciones anteriores, que ya anotaron Park y Kausel [37], tienen gran re-

levancia a la hora de computar numéricamente las integrales. Sugieren en primer lugar


0 5 10 15 20−3

−2

−1

0

1

2

Ω [−]

Kxx

(Ω)s

in(Ω

τ)

[−]

τ = π/2τ = π/2

(a) Bxx

0 5 10 15 20−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

Ω [−]

Kxz

(Ω)s

in(Ω

τ)

[−]

τ=π/2τ=π/2

(b) Bxz

0 5 10 15 20−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

Ω [−]

Kzz

(Ω)s

in(Ω

τ)

[−]

τ=π/2τ=4 π/2

(c) Bzz

Figura 4.4: Evaluación de los núcleos de las integrales Bxx, Bxz, Bzz para el caso concreto de z = 0, ν =0.25.

dividir la integral en dos intervalos: [ΩS,ΩP] y [ΩP,∞], y en segundo lugar aproximarlos núcleos por sus asíntotas respectivas a partir un determinado valor de Ω = Ωa

para el que el error es despreciable. De hecho, Park y Kausel [37] propusieron dividirla integral en tres intervalos [ΩS,ΩP] [ΩP,Ω =

√2] y [ΩP,∞], donde el punto Ω =

√2

coincide con los ceros de los núcleos de Bxx, Bzz, y a partir del cual cambian de signo.Sin embargo, se puede observar que en éste punto no se existe cambio brusco de latendencia de los núcleos, esto es, es derivable y las curvas permanecen suaves. Trasrealizarse varios experimentos numéricos, el hecho de partir adicionalmente las inte-grales no se traduce en un descenso del tiempo de cómputo sino más bien lo contrario,en término medio necesita un 47% más de tiempo. En la Tabla 4.1 se computa el costeadicional en porcentaje al considerar esta nueva partición de las integrales en Ω =

√2

para los tres términos Bxx, Bxz, Bzz y para dos casos correspondientes a ν = 0.25, 0.33.Se ha considerado el receptor sobre la superficie y se ha realizado la media de los resul-tados calculados para τ comprendido entre [π/18, 20π] siendo ∆τ = π/18. Por otrolado, Park y Kausel [37] también propusieron tomar el valor Ωa = 3/a a partir del cualaproximar los núcleos por sus respectivas asíntotas cometiendo un error despreciable.En este sentido en la Tabla 4.2 se muestra el error cometido tras considerar tres valores


0 2 4 6 8 100

0.5

1

1.5

2

2.5

3

Ω [−]

Kxx

(Ω)

[−]

(a) Bxx

0 2 4 6 8 10−0.6

−0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

Ω [−]

Kxz

(Ω)

[−]

(b) Bxz

0 2 4 6 8 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Ω [−]

Kzz

(Ω)

[−]

(c) Bzz

Figura 4.5: Evaluación de los núcleos de las integrales Bxx, Bxz, Bzz para el caso concreto de z = 0, ν =0.25 y obviando el término senoidal (linea continua). En línea discontinua se muestra las asíntotas decada núcleo respectivo.

de Ωa (3/a, 5/a, 7/a) para los términos Bxx, Bxz, Bzz Igualmente se ha considerado doscasos diferentes de ν = 0.25, 0.33, y el receptor sobre la superficie. Se ha realizado lamedia del error para los mismos valores de τ considerados anteriormente. Evidente-mente el error es tanto menor cuanto mayor es el valor de Ωa y, como no puede ser deotra forma, el error para el término Bxz o Bzx es nulo puesto que el núcleo integral seanula para Ωa > 1/a. En este trabajo se ha tomado un valor de aproximación asintó-tica más conservador que Park y Kausel de Ωa = 5/a en el que el error se mantienemenor 0.25%.

Coste adicional computacional (%) Bxx Bxz Bzz

ν = 0.25 54.60 42.66 54.23

ν = 0.33 58.23 18.61 57.81

Tabla 4.1: Coste adicional adicional computacional de la rama integral Bij al considerar tres intervalos deintegración respecto de considerar dos intervalos. C.P.A = 100(C.P.3int. − C.P.2int.)/C.P.2int.


Error (%) Ωa = 3/a Ωa = 5/a Ωa = 7/a

Bij Bxx Bxz Bzz Bxx Bxz Bzz Bxx Bxz Bzz

ν = 0.25 0.36 0 0.15 0.077 0 0.032 0.028 0 0.011

ν = 0.33 1.11 0 6.7× 10−3 0.25 0 4.8× 10−4 0.088 0 3.3× 10−4

Tabla 4.2: Error porcentual de los términos Bxx, Bxz, Bzz al aproximarlos por sus respectivas asíntotas paratres valores diferentes del puto de aproximación Ωa. Se considera el receptor sobre la superficie (z = 0).Los valores están promediados para varios valores de τ comprendidos entre [π/18, 20π].

Tabla 4.3: Funciones asintóticas

Bij =∫ ∞

ΩaKij sinΩτdΩ, Kij ≈ Kij,a

Kij,a z = 0 z 6= 0

Kxx,a − 12π

4Ω

− 1π

(4

Ω4 cos τpΩ + 2Ωcos τsΩ

)

Kxz,a 0 − 1π

(−2Ω2 sin τpΩ + 4

Ω2 sin τsΩ)

Kzx,a 0 − 1π

(4aΩ2 sin τpΩ − 2

Ω2 sin τsΩ)

Kzz,a − 12π

4aΩ

− 1π

(4aΩ3 cos τsΩ + 2a

Ω cos τpΩ)

En suma, la rama integral concerniente a transformada inversa de Fourier (Bij) de-be dividirse al menos en dos intervalos, el primero definido entre los límites [1, 1/a],yel segundo entre [1/a,∞). Además, para éste último intervalo, los núcleos puedenaproximarse por sus asíntotas a partir del valor de aproximación Ωa, con lo cual elintervalo debe subdividirse nuevamente en [1/a,Ωa] y [Ωa,∞). La ventaja de utilizarlas funciones asintóticas es que son muy sencillas y permiten calcular las integrales deforma analítica y cerrada, por lo que se elimina el problema de evaluar numéricamen-te integrales impropias. Estas asíntotas se calculan sin más que tomar el límite de losnúcleos para Ω → ∞. En la Tabla 4.3 se muestran las funciones asintóticas, eliminandoel término senoidal sin τΩ, de cada uno de los núcleos integrales. En cuanto a las inte-grales definidas para los intervalos comprendidos entre [1, 1/a] y [1/a,Ωa], deben serevaluadas numéricamente. Debido a su naturaleza oscilatoria, Park y Kausel proponenutilizar el método de integración Clenshaw-Curtis con polinomios de Chebyshev [49].

Aún con las anotaciones realizadas, la evaluación de la función de Green es costosacomputacionalmente, sobre todo debido a la integración numérica. Sin embargo, gra-cias a la adimensionalización realizada debe observarse que la rama integral Bij(d, τ, ν)solo depende de la profundidad, del tiempo adimensional y del coeficiente de Poissonν. De hecho para el caso particular en el que el punto de evaluacion es superficial, la


dependencia de la profundidad desaparece. Park Kausel [37] propusieron que, paracálculos repetitivos, se podría precomputar las integrales para diferentes valores deν y una densa discretización temporal τ,. De esta forma la evaluación de la rama in-tegral es virtualmente inmediata, e independiente de las propiedades dinámicas delsemiespacio, por rápido acceso a la tabla. Se debe percatar que la rama integral exhibeun periodo principal de oscilación de τ = 2π/3, por lo que debe elegirse un paso tem-poral (∆τ,) suficientemente pequeño para evitar el fenómeno de aliasing. Para valoresintermedios de τ puede realizarse una simple interpolación.En la Figura 4.6 se muestra la evolución temporal de la función de Green P-SV ante

una carga impulsiva tipo senoidal. Se ha considerado un semiespacio con las mismaspropiedades dinámicas que las elegidas en el apartado. La fuente y receptor super-ficiales se separan también 10 m. Se ha utilizado la formulación presentada en esteapartado y se ha realizado la transformada inversa de Fourier sobre el número deonda numéricamente de forma discreta. Además para cada una de las componentes,se ha superpuesto la función de Green en el dominio natural espacio-tiempo exactapropuesta por Eringen y Suhubi [16] para un impulsivo delta de Dirac. Puede apre-ciarse que ambas soluciones prácticamente coinciden, excepto que la singularidad delas ondas de Rayleigh es finita para la solución de Kausel debido a la naturaleza se-noidal finita de la fuente. Recalcar de nuevo que la respuesta se presenta en formaadimensional (GH

ij = Gijµπr/Z) frente al tiempo adimensional τ = tcs/r.

4.3 solución de una placa estratificada mediante el tlm

El "Thin-Layer Method"(TLM), tal como lo introduce su autor Kausel [22], es unatécnica o herramienta semi-analítica que puede ser usada para resolver problemas depropagación de ondas en un medio parcialmente heterogéneo en alguna dirección, porejemplo, una placa estratificada, ante una carga dinámica cualquiera. Aquí se utilizarápara hallar la función de Green ante una carga concentrada. En esencia el método sebasa en una discretización parcial del medio en la dirección estratificada, en donde losdesplazamientos son aproximados por funciones de interpolación, mientras que parael resto de direcciones la solución es exacta . De esta forma se llega a una matriz derigidez que relaciona los desplazamientos con las tensiones aplicadas en el dominionúmero de onda-frecuencia. La respuesta en el tiempo puede computarse numérica-mente mediante una trasformada inversa de Fourier, o mejor, mediante el método dela superposición modal.Históricamente, el método lo introdujo Lysmer [28] en 1970 para estudiar las ondas

de Rayleigh y ondas de Love en un medio horizontalmente estratificado. Sin embargofue Kausel [26, 20] en 1981 quién desarrollo profundamente el método para estudiarproblemas de propagación de ondas en un medio isótropo heterogéneo ante una car-ga dinámica en el dominio transformado número de onda-frecuencia. Posteriormentekausel [21] extendió este método a medios anisótropos, y en 1994 formuló el método

4.3 solución de una placa estratificada mediante el tlm 75

0 0.5 1 1.5 2

−500

0

500

Tiempo [−]

Des

plaz

amie

nto

[−]

(a) Gδxx

0 0.5 1 1.5 2

−500

0

500

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) Gδxz

0 0.5 1 1.5 2

−500

0

500

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) Gδzx

0 0.5 1 1.5 2

−500

0

500

Tiempo [−]

Des

plaz

amie

nto

[−]

(d) Gδzz

Figura 4.6: Evolución temporal de las componentes de la función Green adimensional, Gδ(10, 0, 0, t; 0, 0, 0)en el plano (P-SV) de un semiespacio homogéneo. Linea continua negra corresponde con la función deGreen propuesta por Park y Kausel, y los círculos grises corresponde con la exacta propuesta por Eringeny Suhubi

directamente en el dominio del tiempo para problemas planos ante una carga linealSH y P-SV [22]. Su discípulo, Park [35], en 2002 formuló el método en coordenadascilíndricas (CTLM) para obtener las funciones de Green en un medio estratificado an-isótropo. Y en 2006 Park y Kausel [24, 25] presentaron un método híbrido directamenteel dominio del tiempo para obtener las funciones de Green del semi-espacio infinitohorizontalmente estratificado ante una carga lineal SH y P-SV, y una puntual.

En este apartado se introduce el TLM en el dominio del tiempo y en coordenadascartesianas para obtener la función de Green de una placa estratificada de espesorfinito en la dirección vertical z. Inicialmente, se presenta la base matemática inicial delmétodo [22, 35]. En los siguientes apartados se presenta la formulación particular parael problema en el plano ante una carga lineal P-SV y el antiplano ante carga lineal SHen el dominio del tiempo [24, 25]. Por último ambas soluciones son validadas con lasolución analítica de una placa homogénea presenta por Park [35].


4.3.1 Planteamiento del método en coordenadas cartesianas

Se considera un medio de espesor finito horizontalmente estratificado, anisótropo,lineal y elástico, en el que actúa una carga dinámica. La ecuación de Navier, para unestrato localmente homogéneo, en forma matricial puede expresarse como:

ρu− LTDL = b (4.36)

(4.37)

donde u = [ux, uy, uz]T es el vector desplazamiento, D es la matriz de propiedades ennotación de Voigt , b son las fuerzas de volumen y L es el operador diferencial quepuede escribirse como:

L = Lx∂

∂x+ Ly

∂

∂y+ Lz

∂

∂z(4.38)

donde:

Lx =

1 0 0

0 0 0

0 0 0

0 0 0

0 0 1

0 1 0

,Ly =

0 0 0

0 1 0

0 0 0

0 0 1

0 0 0

1 0 0

,Lz =

0 0 0

0 0 0

0 0 1

0 1 0

1 0 0

0 0 0

(4.39)

Por otro lado, el vector de tensiones internas en un plano horizontal se expresa:

s = LzTDL (4.40)

cuyas componentes son s = [τxz, τyz, σz]T

En cuanto a las condiciones de contorno en tensiones sobre la cara superior y la carainferior se pueden expresar como:

t− s ˚ = 0 (4.41)

siendo t el vector de tracciones externas y s ˚ el vector de tensiones internas en lainterfase de normal ν.Para resolver esta ecuación del movimiento sujeta a las condiciones de contorno en

tensión, Kausel propone dividir el sólido horizontalmente estratificado en delgadas


láminas desde el punto de vista de los elementos finitos, cada una de las cuales sesubdivide a su vez en sub-láminas. Ésta última división atiende al grado m de in-terpolación, de forma que cada lámina (l) se compone de las interfases 1, 2, ...,m+ 1nombradas de arriba a abajo. El equilibro se impone al aplicar tracciones en la super-ficie superior e inferior. En cuanto a los desplazamiento dentro de la lámina, similar alos elementos finitos, se aproxima mediante una función de interpolación en direcciónla dirección de la discretización z.Véase Figura 4.7. Esto es:

Figura 4.7: Discretización de una lámina. Kausel [22].

ul = NUl (4.42)

siendo N(z) una matriz que tiene los polinomios de interpolación, dependiente delgrado de interpolación, y Ul(x, y, t) un vector que contiene los desplazamientos en lasinterfases 1, 2, ...,m+ 1 de la lámina l:

Ul = [ul1T,ul

2T, ...,ul

m+1T]T (4.43)

N = [ξI, (1− ξ)I] para el caso lineal tratado aquí (4.44)

I es la matriz identidad y ξ = z/h la componente vertical adimensional.Introduciendo esta discretización en la Ecuación de Navier (4.36) y en la condición

de contorno (4.41) para una lámina, se llega a un sistema de ecuaciones desequilibra-das en donde, tal como Kausel [22] aclara, es necesario introducir unas fuerzas devolumen residuales rl y de contorno ql para equilibrar dichas ecuaciones:

bl − ρul − LTDL = rl

tl − sl˚ = ql (4.45)

Para obtener la ecuación del movimiento de la lámina, basta con imponer que lasfuerzas residuales producen un trabajo virtual nulo:


[

δul1Tql1 + δul

m+1Tqlm+1 +

∫ h

0δulTrldz

]

= 0 (4.46)

Desarrollando las expresiones y reagrupando términos se llega finalmente a la ex-presión:

pl1

pl2...

plm+1

=∫ h

0NTNρdz

ul1

ul2...

ulm+1

+

LTzDLul|z=h

0...

LTzDLul|z=0

−∫ h

0NTLTDLNρdz

ul1

ul2...

ulm+1

(4.47)

que relaciona las fuerzas exteriores pl en las interfases (engloba las tracciones tl yfuerzas de volumen bl), con las fuerzas internas a través de los desplazamientos en lainterfases ul. De forma matricial y compacta, esta ecuación puede expresarse alterna-tivamente por:

Pl = MlUl −Alxx

∂2Ul

∂x2−Al

xy

∂2Ul

∂x∂y−Al

yy

∂2Ul

∂y2− Bl

x

∂Ul

∂x− Bl

y

∂Ul

∂y+GlUl (4.48)

donde:

Ml =∫ h

0NTNρdz

Alii =

∫ h

0NTDijNρdz

Alxy =

∫ h

0NT(Dxy +Dyx)Nρdz

Bli =

∫ h

0NTDizNρdz−

∫ h

0NTDziNρdz

Gl =∫ h

0N′TDzzN

′ρdz

i = x, y (4.49)

(4.50)

en el que se ha definido:

N′ =∂N

∂z

Dij = LTi DLj, i, j = x, y, z (4.51)


Como puede verse, la Ecuación del movimiento (4.48) se expresa en forma diferen-cial para las coordenadas x e y, mientras que la coordenada z se elimina debido a ladiscretización en esa dirección. Para resolver esta ecuación, Kausel [22] propone rea-lizar la transformada de Fourier sobre el espacio (x, y), suponiendo invarianza en elplano horizontal, y sobre el tiempo t. De esta forma se obtiene una ecuación en formade matriz de rigidez dinámica que relaciona las tracciones en las interfases con losdesplazamientos en el dominio k− ω:

Pl = KlUl (4.52)

siendo

Kl = k2xAlxx + kxkyA

lxy + k2yA

lyy + i(kxB

lx + kyB

ly) +Gl − ω2Ml (4.53)

donde la tilde significa doble trasformación del espacio-tiempo al número de onda-frecuencia.

Hasta aquí se ha presentado la formulación para una sola lámina, donde las ma-triz elemental Kl y por ende las matrices Al,Bl,GlMl, así como los desplazamientosUl y tracciones Pl se refieren a la lámina l. La ecuación del movimiento para un me-dio estratificado se realiza mediante el acoplamiento de las n matriz elementales quecomponen la placa en el sentido de los elementos finitos, esto es:

K = K1 ⊕ K2 . . . Kn−1 ⊕ Kn (4.54)

donde el operador ⊕ significa el acoplamiento de las matrices tal como se representaen la Figura 4.8. Evidentemente este acoplamiento es extrapolable al resto de matriceselementales.

Efectivamente se llega a una ecuación del movimiento que relaciona los desplaza-miento en las N = nm+ 1 interfases con unas tracciones externas genéricas sobre lasmismas para un medio de profundidad finita. En presente trabajo se utiliza este méto-do para obtener la función de Green de una placa isótropa horizontalmente estratifica.En este caso, kausel propuso modificar la componente vertical de los desplazamientoy tracciones por el factor imaginario i, resultado así una matriz de rigidez K real ysimétrica:

P = KU (4.55)

donde ahora el vector desplazamiento, tracciones y matriz de rigidez quedan modifi-cados por:


Figura 4.8: Acoplamiento de n láminas en el sentido de los elementos finitos. Kausel [24].

U = [ux1, uy

1, iuz1, ..., ux

N , uyN , iuz

N ]

P = [px1, py1, i pz1, ..., pxN , pyN , i pzN ]

K = k2xAxx + kxkyAxy + k2yAyy + kxBx + kyBy +G− ω2M (4.56)

y en el que Bx es real y simétrica.

4.3.2 Problema bidimensional anti-plano (SH)

En este apartado se expone la función de Green, obtenida por Kausel [22], de unplaca horizontalmente estratifica en el dominio k − t sujeta a una carga lineal SH.Evidentemente para este problema plano los términos dependientes de la coordenaday deben ser eliminados de la Ecuación (4.55), y debido al desacoplamiento de losdesplazamientos y tracciones solo debe retenerse la componente fuera del plano. Deesta forma la ecuación del movimiento se reduce a:

[k2xAxx +G− ω2M]Uy = Py (4.57)

Las matrices Axx,G y M dependen de las propiedades de los estratos y su estructu-ra depende del grado de interpolación. En la Tabla 4.4 se muestran para una inter-polación lineal, en Kausel [22] y Park [35] pueden econtrarse expresiones para unainterpolación cuadrática. En cuanto a la carga se considera que actúa, sin pérdida degeneralidad, en el origen x = 0, y a una profundidad d que debe coincidir con alguna


Tabla 4.4: Matrices elementales del TLM. Problema plano SH.

Alxx =

µlhl

12

[

5 1

1 5

]

,Gl =µl

hl

[

1 −1

−1 1

]

,Ml =ρlhl

12

[

5 1

1 5

]

de las interfases fruto de la discretización vertical realizada. En el dominio k− ω éstaviene dada por:

Py = δ(zi − d) (4.58)

donde δ(zi − d) es un vector delta de Dirac en la que todas las componentes son ceroexcepto para la profundidad d.

La respuesta en el tiempo se realiza mediante una expansión modal. Los modos depropagación del medio se obtiene sin más que establecer a cero el vector de cargasexternas, y para un determinado número de onda, resolver un problema estándarlineal de autovalores y autovectores:

[k2xAxx +G︸︷︷︸

Kestatica

−ω2jM]φj = 0 (4.59)

que en forma matricial se expresa:

[Kestatica −MΩΦj] = 0 (4.60)

donde ωj representa los autovalores y φj los autovectores asociados o modos de propa-gación. En forma matricial se expresan como Ω = diag[ω] y Φ = [φj]. Adicionalmente,como aclara Kausel [22], los autovectores deben de satisfacer las siguientes condicionesde ortogo-normalidad:

ΦMΦT = 0

ΦKestaticaΦT = −Ω

2 (4.61)

Notar que la matrices Axx y M son definida positiva, mientras que G es semidefi-nidas positivas. Por tanto todos los autovalores deben ser reales y no negativos. Dehecho para el caso particular de Kx = 0 existe un movimiento como solido rígido conω0 = 0 y todas las componentes del autovector asociado φ0 = cte..Para obtener la solución modal se supone que la respuesta es de la forma:


Uy = Φf (4.62)

Intoduciéndolo en la Ecuación del movimiento (4.57) se tiene que la función f es iguala:

f = [Ω2 − Iω2]−1Φ

Tδ(zi − d) (4.63)

Sustituyendo en la Ecuanción (4.62), los desplazamiento en el dominio k−w se expresacomo:

Uy = Φ[Ω2 − Iω2]−1Φ

Tδ(zi − d) (4.64)

Finalmente, la respuesta en el dominio k− t se obtiene realizando una transformadainversa de Fourier sobre el número de onda:

Gδy(kx, t) = Uy(kx, t) =

∫ ∞

−∞Uy(kx, t)eiωt = ΦHΦ

Tδ(zi − d) (4.65)

siendo H = diag[sin ωjt/ωj] y el símbolo ˇ sobre la variable representa la trasformadasimple al número de onda.Alternativamente,tal y como presenta Kausel en su paper [24], esta ecuación puede

expresarse en forma de sumatorio:

Gδ,edy =

N

∑j=1

sinωjt

ωjφj

eφjd (4.66)

donde e se refiere al índice de la interfase a la profundidad del receptor, y d a la de lafuente.Notar que la solución es valida para un medio sin atenuación. Igual que para el caso

del semiespacio presentado en el apartado 4.2, es posible introducirle amortiguamien-to sin más que realizarle la siguiente modificación:

Gδ,edy =

N

∑j=1

e−ξ jωjsinωdjt

ωdjφj

eφjd (4.67)

donde ξ j representa el factor de amortiguamiento critico modal, que puede depender

del número de onda y ωdj = ωj

√

1− ξ2j


En la Figura 4.9 se muestra una comparativa de la solución obtenida con el TLMde una placa homogénea de espesor 5 m respecto de la solución analítica. Las pro-piedades del medio son: ρ = 2000kg/m3 , cs = 250 m/s, cp = 467.70 m/s . Se haconsiderado una carga impulsiva tipo senoidal en dirección y superficial y se ha dibu-jado los desplazamientos fuera del plano para en un punto situado a 4 m de distanciade la fuente sobre la superficie. Los desplazamientos se han adimensionalizados poruyyµxπ/X siendo X la magnitud de la fuerza, mientras que el tiempo por t/(x/cs).La placa se ha discretizado en 140 elementos y se han tomado los 140 primeros modosde vibración. En la gráfica puede verse el tiempo de llega de las ondas S para tS = 1,tras el cual se registra la reflexión de las ondas en la cara inferior del plano tSS = 2.7.Además puede percatarse que ambas curvas son prácticamente idénticas.

0 1 2 3 4 5−200

0

200

400

600

800

1000

1200

1400

1600

Tiempo [−]

Des

plaz

amie

nto

[−]

Figura 4.9: Evolución temporal de las componentes de la función Green adimensional, Gyy(5, 0, 0, t; 0, 0, 0)

en el plano (SH) de una placa homogénea. En línea continua se dibuja la solución del TLM y en círculosla solución exacta

4.3.3 Problema bidimensional en el plano (P-SV)

Continuando con el método el TLM, a continuación se presenta la función de Greende una placa horizontalmente estratifica ante una carga lineal P-SV [25]. La ecuacióndel movimiento es análoga al caso anterior SH, excepto que ahora se retienen lascomponentes en el plano de los desplazamientos y las tracciones, lo que significa queahora éstas sean vectores de dos componentes relativos a la dirección horizontal yvertical. Además aparece una nueva matriz Bx que tiene en cuenta el acoplamiento deambas direcciones. La ecuación del movimiento se reduce así pues a:

[k2xAxx + kxBx +G− ω2M]U = P (4.68)


Tabla 4.5: Matrices elementales del TLM. Problema plano P-SV.

Alxx =

hl

12

5(λl + 2µl) 0 (λl + 2µl) 0

0 5µl 0 µl

(λl + 2µl) 0 5(λl + 2µl) 0

0 µl 0 5µl

,Blx =

12

0 −(λl − µl) 0 (λl + µl)

−(λl − µl) 0 −(λl + µl) 0

0 −(λl + µl) 0 (λl − µl)

(λl + µl) 0 (λl − µl) 0

Gl = 1hl

µl 0 −µl 0

0 λl + 2µl 0 −(λl + 2µl)

−µl 0 µl 0

0 −(λl + 2µl) 0 λl + 2µl

,Ml =ρlhl

12

5 0 1 0

0 5 0 1

1 0 5 0

0 1 0 1

donde:

U = [Ux, iUz]

P = [Px, iPz] (4.69)

En la Tabla 4.5 se muestran el contenido de las matrices elementales para una inter-polación lineal, en Park [35] pueden encontrarse expresiones para una interpolacióncuadrática. En cuanto a la carga se considera que se aplica en el origen x = 0 y a unaprofundidad d, e indistintamente puede actuar en la dirección horizontal o vertical.Además se multiplica por la matriz Q que convierte la parte imaginaria en puramentereal, lo que significa que los resultados sean también puramente reales. En el dominiok− ω ésta viene dada por:

Pℜ = QP, Q =

[

1 0

0 −i

]

(4.70)

Pℜ = [δ(zi − d), 0]T dirección horizontal

Pℜ = [0, δ(zi − d), ]T dirección vertical (4.71)

(4.72)

donde δ(zi − d) es un vector delta de Dirac de longitud N = nm+ 1 en la que todaslas componentes son cero excepto para la profundidad d, y 0 representa un vector deceros.La respuesta en el tiempo nuevamente se realiza mediante una superposición modal.

Los modos de propagación del medio se obtiene al establecer a cero el vector de cargasexternas, y para un determinado número de onda, se resuelve un problema estándarlineal de autovalores y autovectores:

[k2xAxx + kxBx +G︸︷︷︸

Kestatica

−ω2jM]φj = 0 (4.73)


que en forma matricial se expresa:

[Kestatica −MΩΦj] = 0 (4.74)

donde ωj representa los autovalores y φj los autovectores asociados o modos de pro-pagación. En forma matricial se expresan como Ω = diag[ω] y Φ = [φj]. Fíjese queen este caso deben de existir el doble de autovalores y autovectores respecto al pro-blema fuera del plano j = 1, 2...2N debido a que se consideran los desplazamientoen las dos direcciones del plano. Nuevamente los autovectores deben de satisfacer lascondiciones de ortogo-normalidad, Kausel [22]:

ΦMΦT = 0

ΦKestaticaΦT = −Ω

2 (4.75)

Igualmente, se deduce que la matriz de rigidez es definida o semidefinida positivaen función del valor del número de onda, por lo que todos los autovalores deben serreales y no negativos. Además para caso particular de Kx = 0 se recogen dos movi-mientos de solido rígido asociado a las dos direcciones en el plano, esto es, ω0,ω1 = 0.

Nuevamente, para obtener la solución modal se supone que la respuesta es de laforma:

Uℜ = Φf (4.76)

Intoduciéndolo en la Ecuación del movimiento (4.68) se tiene que la función f es iguala:

f = [Ω2 − Iω2]−1Φ

TPℜ (4.77)

Sustituyendo en la Ecuanción (4.76), los desplazamiento en el dominio k−w se expresacomo:

Uℜ = Φ[Ω2 − Iω2]−1Φ

TPℜ (4.78)

Igualmente, realizando la transformada inversa de Fourier sobre el número de onda,la respuesta en el dominio k− t es:

Uℜ(kx, t) =∫ ∞

−∞Uℜ(kx, t)eiωt = ΦHΦ

TPℜ (4.79)


siendo H = diag[sin ωjt/ωj] y el símbolo ˇ sobre la variable representa la trasformadasimple al número de onda..En este punto, se debe notar que los autovalores y autovectores pueden dividirse en

dos relativos a las dos direcciones en el plano, y las tracciones pueden actuar indepen-dientemente también en cualquier dirección. Por tanto la función de Green puramentereal en el dominio k− t en forma cerrada puede expresarse como:

Gδℜ(k, t) =

[

ΦxHΦxT[δ(zi − d), 0]T ΦxHΦ

Tx [0, δ(zi − d)]T

ΦzHΦzT[δ(zi − d), 0]T ΦzHΦ

Tz [0, δ(zi − d)]T

]

(4.80)

Nuevamente,a esta función de Green se le puede introducir un amortiguamiento enel mismo sentido que para la función de Green SH. Así pues en forma de sumatorio,esta se expresa como:

Gδ,edij,ℜ =

N

∑j=1

e−ξ jωjsinωdjt

ωdj

[

φxjeφxj

d φxjeφzj

d

φzjeφxj

d φzjeφzj

d

]

(4.81)

en el que se ha definido:

ωdj = ωj

√

1− ξ2j (4.82)

donde e y d representan la profundidad del receptor y de la fuente respectivamente,y ξ j representa el factor de amortiguamiento critico modal, que puede depender delnúmero de onda.Por último La función de Green actual se obtiene sin más que revertir la trasforma-

ción imaginaria realizada anteriormente, esto es:

Gδ(k, t) = Q−1Gδℜ(k, t)Q (4.83)

En la Figura 4.10 se muestra una comparativa de la solución obtenida con el TLM deuna placa homogénea de espesor 5 m respecto de la solución analítica. Las propieda-des del medio son: ρ = 2000kg/m3 , cs = 250 m/s, cp = 467.70 m/s . Se ha consideradouna carga impulsiva tipo senoidal vertical superficial y se ha dibujado los desplaza-mientos verticales para en un punto situado a 4 m de distancia de la fuente sobre lasuperficie. Los desplazamientos se han adimensionalizados por uzzµxπ/X siendo X

la magnitud de la fuerza, mientras que el tiempo por t/(x/cs). La placa se ha discre-tizado en 140 elementos y se han tomado los 140 primeros modos de vibración. En lagráfica puede observase la llega de las ondas P , S y Rayleigh, tras los cuales apareceuna perturbación debido a las ondas reflejadas en la interfase inferior. Además puedepercatarse que ambas curvas son prácticamente idénticas.

4.4 solución del semiespacio estratificado 87

0 0.5 1 1.5 2−1000

−500

0

500

1000

1500

Tiempo [−]

Des

plaz

amie

nto

[−]

Figura 4.10: Evolución temporal de las componentes de la función Green adimensional,Gzz(5, 0, 0, t; 0, 0, 0) en el plano (P-SV) de una placa homogénea. En línea continua se dibuja la solucióndel TLM y en círculos la solución exacta

4.4 solución del semiespacio estratificado

En las secciones precedentes se ha presentado la función de Green de una placahorizontalmente estratifica, y la función de Green del semiespacio homogéneo parael problema plano. Siguiendo, la guía introducida en este capítulo, el siguiente pasoconsiste en acoplarlas para obtener la función de Green del semiespacio estratificado.Este acoplamiento propuesto por Park y Kausel [24, 25], recibe el nombre de métodohíbrido. En esta sección, en primer lugar se presenta el método híbrido para el proble-ma plano ante una carga lineal SH y P-SV. Posteriormente se extrapola al problematridimensional ante una carga puntual.

4.4.1 Problema bidimensional

En primer lugar es conveniente aclarar la nomenclatura de la que se va a hacer uso.Aunque para el problema plano, los desplazamientos pueden desacoplarse, en estasección se considera el tensor de Green completo (3× 3) tanto para la placa comopara el semiespacio. Además para hacer distinción de cada una de ellas, en adelante elsuperíndice SE hace referencia al semiespacio y el superíndice P a la placa, y se omiteel superíndice δ que se refiere al tipo de carga delta de Dirac. Matemáticamente ambassoluciones se expresan como siguen:


GSE(kx, z, t, z′) =

GSExx 0 GSE

xz

0 GSEyy 0

GSEzx 0 GSE

zz

función de Green del semiespacio (4.84)

GP(kx, z, t, z′) =

GPxx 0 GP

xz

0 GPyy 0

GPzx 0 GP

zz

función de Green de la placa (4.85)

Se considera un semiespacio compuesto por una serie de estratos horizontales bajolos cuales se extiende un medio homogéneo e infinito. La idea del método consisteen dividir el semiespacio en dos regiones, una región finita que engloba a todos losestratos, y otra región infinita homogénea. Ambas regiones se tratan de forma inde-pendiente tal que la primera se modela mediante la solución de la placa y la segundamediante la del semiespacio homogéneo. El método se formula de forma diferentesegún se encuentre la carga en la región estratifica o en el semiespacio, aunque eldesarrollo es análogo. Se supone el caso que la carga está dentro de la placa, y se dejaal lector su extrapolación al otro caso. Entonces se impone a la placa el equilibrio, loque requiere que en la superficie inferior existan unas tensiones τ(kx, t) en principiodesconocidas. Evidentemente estas tensiones son debidas a las reacciones del semies-pacio sobre los estratos, tanto que en la superficie superior del semiespacio actúan lasmismas tensiones pero de signo contrario. El acoplamiento entre la placa estratifica yel semiespacio se realiza sin más que aplicar el principio de superposición, lo que per-mite calcular las tracciones en la interfase τ(kx, t). Conocidas éstas es sencillo calcularlos desplazamientos en cualquier región. Matemáticamente, paso por paso, el métodose expresa tal como sigue. Se ha tomado como referencia del sistema de coordenadasla interfase de unión, donde el vector ~ez apunta hacia fuera del semiespacio (véaseFigura 4.11):

Figura 4.11: Superposición de la respuesta de la placa estratifica con el semiespacio homogéneo. Kausel[24].


El campo de desplazamiento en cualquier lugar de la placa estratifica viene dadopor el producto de convolución del tensor de Green por la carga externa más elproducto de convolución del tensor con las tensiones en la interfase, a prioridesconocidas:

U(kx, z, t, z′) = GP(kx, z, t, z′) ∗ P(kx, t)− GP(kx, z, t, 0) ∗ τ(kx, t) (4.86)

Que particularizando para la cara inferior se expresan:

U(kx, 0, t, z′) = GP(kx, 0, t, z′) ∗ P(kx, t)− GP(kx, 0, t, 0) ∗ τ(kx, t) (4.87)

Análogamente, el campo de desplazamiento en el cualquier lugar del semiespa-cio viene dado por:

U(kx, z, t, 0) = GSE(kx, z, t, 0) ∗ τ(kx, t) (4.88)

Y en la interfase:

U(kx, 0, t, 0) = GSE(kx, 0, t, 0) ∗ τ(kx, t) (4.89)

Evidentemente, los desplazamiento en la interfase común tienen que ser igualespara ambas regiones, con lo que aplicando compatibilidad se llega al siguientesistema de ecuaciones:

GP(kx, 0, t, z′) ∗ P(kx, t) = [GSE(kx, 0, t, 0) + GP(kx, 0, t, 0)] ∗ τ(kx, t) (4.90)

Que en la forma discreta de sumatorio se puede expresar como:

i

∑j=0

Hi−jPj =i

∑j=0

Fi−jτ j (4.91)

donde se han definido las matrices:

H = GP(kx, 0, t, z′) (4.92)

F = GSE(kx, 0, j∆t, 0) + GP(kx, 0, j∆t, 0) (4.93)


Despejando las tracciones τ(kx, t) se llega al siguiente sistema recursivo, en elque debe resolverse secuencialmente debido a la interdependencia con las trac-ciones en los instantes previos:

τ0 = F−10 [H0P0]

τ1 = F−10 [H1P0 +H0P1 − F1τ0]

...

τn = F−10 [HnP0 +Hn−1P1 + ...+H0Pn − (FnP0 + Fn−1P1 + ...+ F1Pn−1)]

(4.94)

Por tanto una vez obtenido las tensiones internas, para obtener el campo dedesplazamiento en cualquier lugar basta con sustituir τ(kx, t) en la Ecuación(4.86) o en la Ecuación (4.88) según la posición del receptor se encuentre dentrode los estratos o en el semiespacio respectivamente.

Evidentemente la función de Green del semiespacio plano estratifico en el dominiok− t se obtiene sin más que aplicar una carga impulsiva en las tres direcciones coor-denadas. Para obtener la respuesta en el dominio espacio-tiempo es necesario realizaruna trasformada inversa de Fourier sobre el número de onda. La ventaja del méto-do precisamente es que, a diferencia del método de la Matriz de Rigidez [26], éstaantitransformada carece de singularidades.

GSEE(kx, z, t, z′) =

δ(t)δ(z−z′)~ex︷︸︸︷

GSEExx

δ(t)δ(z−z′)~ey︷︸︸︷

0

δ(t)δ(z−z′)~ez︷︸︸︷

GSEExz

0 GSEEyy 0

GSEEzx 0 GSEE

zz

(4.95)

en donde, efectivamente la respuesta en el plano está desacoplada de la respuestafuera de él.

GSEE(x, z, t, z′) =∫ ∞

−∞GSEE(kx, z, t, z′)e−ikxxdkx (4.96)

4.4.2 Problema tridimensional

En la sección anterior se ha obtenido la función de Green bidimensional de un se-miespacio horizontalmente estratifica ante una carga lineal SH y P-SV. Como se haintroducido en el apartado 2, la solución tridimensional puede formarse medianteuna combinación de problemas planos. Partiendo de la solución tridimensional en


coordenadas cilíndricas, si se realiza una expansión de Fourier en la coordenada azi-mutal, seguido de una trasformada de Hankel sobre la coordenada radial se obtienela solución en el dominio del numero de onda radial-tiempo. Y efectivamente resultaque ésta solución es virtualmente análoga a la combinación en el número de onda delproblema plano SH y P-SV. Por tanto, combinandado ambas soluciones y realizandouna trasformada inversa sobre el numero de onda radial, se puede calcular la solucióntridimensional en el dominio espacio-tiempo [24, 25, 26].

Supóngase un campo de desplazamiento (o tracciones) tridimensional u(r, θ, t) encoordenadas cilíndricas, la transformación al número de onda y viceversa matemática-mente se expresa así:

u(kr , t) = an

∫ ∞

0rCn

∫ 2π

0Tnu(r, θ, t)dθdr (4.97)

u(r, θ, t) =∞

∑n=0

Tn

∫ ∞

0krCnu(kr , t)dkr (4.98)

en el que se han definido:

an =

1/(2π) n = 0

1/(2π) n = 0(4.99)

Tn =

diag[cos nθ,− sin nθ, cos nθ] caso simétrico

diag[sin nθ, cos nθ, sin nθ] caso antisimétrico(4.100)

Cn =

J′n(krr)nkrr

J′n(krr) 0nkrr

J′n(krr) J′n(krr) 0

0 0 −Jn(krr)

, J′n =

dJn(υ)

υ, υ = krr (4.101)

donde Jn(υ) es la función de Bessel de orden n que satisface las siguientes relacionesde recurrencia:

J − n− 1(krr) + Jn+1(krr) =2nkrr

Jn(krr) (4.102)

J − n− 1(krr)− Jn+1(krr) = 2J′n(krr) (4.103)

Para abreviar, en adelate se denotará Hn a la trasformación de Hankel unidimensio-nal y H−1

n a la trasformación inversa. Matemáticamente aplicado sobre una funcióngenérica representa:

Hn[ f (r)] =∫ ∞

0rJn(krr) f (r)dr (4.104)

H−1n [ f (kr)] =

∫ ∞

0kr Jn(krr) f (kr)dk (4.105)


4.4.2.1 carga en la dirección x

Se considera un semiespacio tridimensional estratificado sujeto a una carga P(x, y, z, t, z′)~exen dirección x impulsiva concentrada aplicada a la profundidad z′, y en el que se co-noce la función de Green bidimensional GSEE

ij (kx, t), que para abreviar en adelante seomite el superíndice.La carga en el dominio espacio tiempo en coordenadas cartesianas viene dada por:

P(x, y, z, t, z′) =

1

0

0

δ(x)δ(y)δ(z− z′)δ(t) (4.106)

y en coordenadas cilíndricas:

P(r, θ, z, t, z′) =

1

0

0

δ(r)δ(θ)δ(t)

r(4.107)

Si este vector es transformado al número de onda-tiempo aplicando la transforma-ción antes detallada (Ecuación (4.97)) y usando la definición simétrica de Tn, se tiene:

P(kr , n, z, t, z′) = an

∫ ∞

0rCn

∫ 2π

0TnP(r, θ, z, t, z′)dθdr =

12π δ1nδ(z− z′)12π δ1nδ(z− z′)

0

(4.108)

donde δ1n es la delta de Kronecker que solo es distinto de cero para n = 1.Es evidente que la carga puntual tridimensional es una combinación de la carga

lineal SH y P-SV en el dominio número de onda-tiempo, en donde ahora el númerode onda se refiere en la dirección radial. En consecuencia, la función de Green delsemiespacio estratifica tridimensional en coordenadas cilíndricas en el dominio (k− t),responde a la misma combinación de las funciones de Green SH y P-SV. Esto es:

Green tridimensional︷︸︸︷

Grx(kr , z, t, z′)

Gθx(kr, z, t, z′)

Gzx(kr, z, t, z′)

=12π

Green bidimensional︷︸︸︷

Gxx(kr , z, t, z′)

Gyy(kr, z, t, z′)

Gzx(kr, z, t, z′)

(4.109)

La función de Green tridimensional del semiespacio estratificado en el dominio delespacio y tiempo se obtiene al aplicar la transformada inversa de Hankel sobre elnúmero de onda con la definición simétrica de Tn [25]:


Grx(r, θ, z, t, z′) =12π

[12H−1

0 [Gxx + Gyy]−12H−1

2 [Gxx − Gyy]] cos θ (4.110)

Gθx(r, θ, z, t, z′) = − 12π

[12H−1

0 [Gxx + Gyy] +12H−1

2 [Gxx − Gyy]] sin θ (4.111)

Gzx(r, θ, z, t, z′) = − 12π

H−11 [iGzx] cos θ (4.112)

(4.113)

4.4.2.2 carga en la dirección y

Ahora se considera la misma carga P(x, y, z, t, z′)~ey en dirección y sobre el semies-pacio estratificado, y donde se conoce igualmente el tensor de Green bidimensionalGij(kx, t).La carga en el dominio espacio tiempo en coordenadas cartesianas se expresa:

P(x, y, z, t, z′) =

0

1

0

δ(x)δ(y)δ(z− z′)δ(t) (4.114)


P(r, θ, z, t, z′) =

0

1

0

δ(r)δ(θ)δ(t)

r(4.115)

Nuevamente si se transforma al número de onda-tiempo pero usando en este casola definición antisimétrica de Tn, se tiene:

P(kr, n, z, t, z′) = an

∫ ∞

0rCn

∫ 2π


12π δ1nδ(z− z′)12π δ1nδ(z− z′)

0

(4.116)

donde δ1n es la delta de Kronecker que solo es distinto de cero para n = 1.Igualmente se deduce que la carga puntual tridimensional es una combinación de

la carga lineal SH y P-SV en el dominio número de onda-tiempo, tal que la función deGreen del semiespacio estratifica tridimensional en coordenadas cilíndricas respondea la misma combinación, esto es:



Gry(kr , z, t, z′)

Gθy(kr, z, t, z′)

Gzy(kr, z, t, z′)

=12π


Gxx(kr , z, t, z′)

Gyy(kr, z, t, z′)

Gzx(kr, z, t, z′)

(4.117)

La función de Green tridimensional del semiespacio estratificado en el dominio delespacio y tiempo se obtiene al aplicar la transformada inversa de Hankel sobre elnúmero de onda con la definición antisimétrica de Tn [25]:

Gry(r, θ, z, t, z′) =12π

[12H−1

0 [Gxx + Gyy]−12H−1

2 [Gxx − Gyy]] sin θ (4.118)

Gθy(r, θ, z, t, z′) =12π

[12H−1

0 [Gxx + Gyy] +12H−1

2 [Gxx − Gyy]] cos θ (4.119)

Gzy(r, θ, z, t, z′) = − 12π

H−11 [iGzx] cos θ (4.120)

(4.121)

4.4.2.3 carga en la dirección z

Por último, se considera la misma carga P(x, y, z, t, z′)~ez en dirección z sobre el semi-espacio tridimensional estratificado, y donde se conoce igualmente el tensor de Greenbidimensional Gij(kx, t).La carga en el dominio espacio tiempo en coordenadas cartesianas se expresa:

P(x, y, z, t, z′) =

0

0

1

δ(x)δ(y)δ(z− z′)δ(t) (4.122)


P(r, θ, z, t, z′) =

0

0

1

δ(r)δ(θ)δ(t)

r(4.123)

Nuevamente, si se transforma al número de onda-tiempo usando la definición simé-trica de Tn, se tiene:

P(kr , n, z, t, z′) = an

∫ ∞

0rCn

∫ 2π


0

0

− 12π δ0nδ(z− z′)

(4.124)

4.5 precisión y convergencia de la solución 95

donde δ0n es la delta de Kronecker que solo es distinto de cero para n = 0.En este caso la carga puntual tridimensional vertical coincide con la carga lineal ver-

tical P-SV multiplicada por el factor − 12π δ0n en el dominio número de onda-tiempo.

Por tanto la función de Green del semiespacio estratifica tridimensional en coorde-nadas cilíndricas proviene de la misma modificación realizada a la función de Greenbidimensional, esto es:


Grz(kr , z, t, z′)

0

Gzz(kr, z, t, z′)

=12π


−Gxz(kr , z, t, z′)

0

−Gzz(kr, z, t, z′)

(4.125)

La función de Green tridimensional del semiespacio estratificado en el dominio delespacio y tiempo se obtiene al aplicar la transformada inversa de Hankel sobre elnúmero de onda con la definición simétrica de Tn [25]:

Grz(r, θ, z, t, z′) = − 12π

H−10 [iGxz] (4.126)

Gθz(r, θ, z, t, z′) = 0 (4.127)

Gzz(r, θ, z, t, z′) =12π

H−10 [iGzz] (4.128)

(4.129)

4.5 precisión y convergencia de la solución

Tal como el propio Kausel [22] describe, el método híbrido es un método semi-analítico en el sentido que se basa en una discretización espacial del semiespacio en ladirección heterogénea. Además, la función de Green en el dominio temporal se obtienemediante una superposición modal, que teóricamente debiera ser infinita, pero en lapráctica tiene que ser truncada para un cierto valor máximo de la frecuencia natural.Y por otro lado, el acoplamiento entre los estratos y el semiespacio se realiza de formadiscreta, lo que requiere elegir un cierto paso temporal.

En éste apartado se estudia la estabilidad y la precisión del método. La estabilidad,como puede intuirse, depende casi exclusivamente del paso temporal elegido para re-solver el sistema de ecuaciones que acopla el semiespacio con la placa estratificada. Encuanto a la precisión, el TLM es el principal responsable, y depende de la discretiza-ción vertical realizada. Éste método es tanto más preciso cuanto menor es el espesorde los elementos tomado evidentemente. De hecho los modos y frecuencias naturalescalculados mediante el método no son exactos, sino presentan una dispersión que seacentúa para frecuencias altas, y depende precisamente de la discretización vertical.


En este sentido Park [] introdujo el parámetro Nλ que significa el número de elemen-tos por longitud de onda, e informa sobre la precisión de la discretización elegida enrelación con la frecuencia natural máxima.A continuación se procede a discutir sobre la estabilidad y precisión en base a una

serie de ejemplos numéricos en los que se juega principalmente con el paso temporaly el número de elementos por longitud de onda. Para ello se considera un semiespaciohomogéneo cuyas propiedades dinámicas son: ρ = 1500kg/m3 , cs = 1000 m/s, cp =

2.08 m/s. A efectos de probar el método se considera virtualmente un estrato de 1000mde profundidad. Además se supone una carga vertical superficial de la forma de unimpulso senoidal tanto temporal como espacial, esto es:

P(r, θ, z, t, z′) = f (r)h(t)δ(z)~ez (4.130)

donde:

f (r) =1

πa2e−(r/a)2 ∀ r ≥ 0 (4.131)

h(t) =2tp

sin2 πt

tp∀0 ≤ t ≤ td (4.132)

donde se ha tomado los periodos tp = 0.1 seg. y a = 100mRealizando la trasformada de Fourier a la función temporal, y la transformada de

Hankel a la función espacial se tiene:

f (kr) =12π

e−(kra/2)2 (4.133)

h(ω) =sinωtp/2

ωtp/2(π2 − (ωtp/2)2)(4.134)

Puede observarse que la amplitud de la función h(ω) decrece con el cubo de lafrecuencia, tal que para ωmax ≥ (4π/tp) es prácticamente despreciable. E igualmentesucede con la función f (kr) que decae exponencialmente y siendo despreciable parakr,max ≥ (2π/a).Parece razonable que para decidir sobre el paso temporal idóneo, se estudie sepa-

radamente los requisitos o condicionantes temporales que ambas soluciones que con-forman el método requieren. Evidentemente la frecuencia modal máxima de corte delTLM debe superior a la frecuencia máxima de la fuente fM ≥ fmax, y el paso temporalelegido no puede ser menor que la frecuencia de Nyquist. Matemáticamente esto seexpresa como:

fmax =2π

ωmax≤ fM ≤ fNyq =

12∆t

(4.135)


Por otro lado, como se comentó en el apartado 4.2.2, el semiespacio homogéneoexhibe una cadencia periódica con el tiempo adimensional τ = 2π

3 , tanto que el pasotemporal debe estar en consonancia con ello:

τSE =2π

3= Tcskmax ⇒ fSE =

3cskmax

2π≤ fNyq =

12∆t

(4.136)

En la Figura 4.12 se muestran los desplazamientos verticales superficiales para tresdiferentes pasos de tiempos: ∆t = 1/(2 fmax), ∆t = 1/(4 fmax) y ∆t = 1/(8 fmax). Enlínea continua se muestra el método híbrido, y en círculos la respuesta exacta delsemiespacio. Los desplazamientos se encuentran adimensionalizados por uzzµrπ/Xsiendo r el radio de separación y X la magnitud de la fuerza, mientras que el tiempopor t/(t/cs). Tal como Park [35] hace notar, el método se vuelve inestable para un pasotemporal ∆t = 1/(2 fmax). De hecho, éste paso es mayor que el que necesita el propiosemiespacio (Ecuación 4.136), y de ahí su inestabilidad. Para el paso temporal inme-diatamente más pequeño ∆t = 1/(4 fmax) parece que el método es estable pero aúnasí parece encontrarse en el límite puesto que coincide exactamente con la frecuenciade Nyquist del semiespacio. En consecuencia, Kausel [25] recomienda tomar un pasotemporal igual a:

∆ ≤ min

[1

8 fmax,

2π

12cskr,max

]

(4.137)

En cuanto a la precisión del método híbrido, depende en mayor medida de la exac-titud de los modos y frecuencias naturales computados mediante el TLM. Es evidenteque la exactitud es tanto mayor cuanto mejor es la discrectización vertical, y obviamen-te cuanto mayor sea también el grado de interpolación utilizado. En este trabajo se hautilizado una interpolación lineal luego no es un parámetro a estudiar. Park [35] en suPhD introdujo con buen criterio la el parámetro Nlambda = λmin/h que significa el nú-mero de elementos o discretizaciones por longitud de onda, siendo λmin = 2π/wmax

y h el espesor de los elementos. Obviamente cuanto mayor sea Nlambda con mayorexactitud se se registra el movimiento de propagación, y en el mismo sentido de lafrecuencia de Nyquist debe ser al menos superior a 2, esto es:

1λmin

≤ fNyq =12h

⇒ Nλ =λmin

h≥ 2 (4.138)

En la Figura 4.13 se muestran los desplazamientos verticales adimensionales paratres discretizaciones diferentes relativas a: Nλ = 2,Nλ = 4 y Nλ = 8. Nuevamenteen línea continua se dibuja el método híbrido y en círculos grises el desplazamientoexacto. Efectivamente, para la peor discretización se aprecia una pérdida de precisiónque se hace más evidente con la llegada de las ondas de Rayleigh. Para la discretizaciónsiguiente se puede observar que prácticamente coincide con la solución exacta y para


0 0.5 1 1.5 2

−2

−1

0

1

2

3

4

Tiempo [−]

Des

plaz

amie

nto

[−]

(a) ∆t = 1/(2 fmax)

0 0.5 1 1.5 2

−2

−1

0

1

2

3

4

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) ∆t = 1/(4 fmax)

0 0.5 1 1.5 2

−2

−1

0

1

2

3

4

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) ∆t = 1/(8 fmax)

Figura 4.12: Estabilidad del método híbrido para tres pasos temporales ∆t = 1/(2 fmax), ∆t = 1/(4 fmax)y ∆t = 1/(8 fmax). Función de Green del semiespacio ante una carga vertical superficial tipo impulsivasenoidal. Gzz(1000, 0, 0, t, 0, 0, 0). En línea continua se dibuja la respuesta del método híbrido y en círculosgrises la respuesta exacta.

Nλ = 8 el error es despreciable. Nótese que en las dos primeras discretizaciones Nλ =

2,Nλ = 4, para un tiempo de t = 1.65 las curvas presentan una pequeña oscilaciónque no debiera aparecer,y se debe presumiblemente al acoplamiento del estrato con elsemiespacio. Puede verse que este efecto desaparece al elegirse una discretización másfina Nλ = 8. Tras estas observaciones, Park [35] recomienda tomar Nλ ≥ 8.


0 0.5 1 1.5 2

−2

−1

0

1

2

3

4

Tiempo [−]

Des

plaz

amie

nto

[−]

(a) Nλ = 2

0 0.5 1 1.5 2

−2

−1

0

1

2

3

4

Tiempo [−]

Des

plaz

amie

nto

[−]

(b) Nλ = 4

0 0.5 1 1.5 2

−2

−1

0

1

2

3

4

Tiempo [−]

Des

plaz

amie

nto

[−]

(c) Nλ = 8

Figura 4.13: Precisión del método híbrido para tres discretizaciones Nλ = 2,Nλ = 4 y Nλ = 8. Función deGreen del semiespacio ante una carga vertical superficial tipo impulsiva senoidal. Gzz(1000, 0, 0, t, 0, 0, 0).En línea continua se dibuja la respuesta del método híbrido y en círculos grises la respuesta exacta.

5E JEMPLOS NUMÉR ICAS

5.1 introducción

Una vez mostrada la formulación de la solución fundamental del semiespacio ho-mogéneo y horizontalmente estratifico, en esta sección se presentan algunos ejemplosnuméricos de interés ingenieril de los que pueden extraerse algunas conclusiones in-teresante. En primer lugar, se cree conveniente validar la implementación de la formu-lación de la solución del semiespacio estratificado con resultados experimentales. Paraello se cuenta con los resultados de un estudio vibratorio del tramo Palencia-León dela futura LAV Madrid-Asturias. Se comparan las curvas de dispersión de las ondassuperficiales experimentales obtenidas en los ensayos de caracterización dinámico delsuelo con las calculadas teóricamente mediante un modelo numérico. Seguidamentese realiza una comparativa entre la función de Green del semiespacio y la del espaciocompleto. Concretamente se estudia la influencia que ejerce la superficie libre del se-miespacio a medida que tanto la fuente como el receptor aumentan de profundidad.En otro apartado se estudia el efecto de introducir amortiguamiento en la respuestadel semiespacio. Y en los dos últimos apartados, se estudia la propagación de ondasen un suelo estratificado para dos casos particulares: en uno se considera un sueloformado por un estrato blando, que simula un terreno vegetal, sobre un semiespaciorocoso; y en el otro se considera un suelo blando en el que se le intercala un fino es-trato rocoso. Se analiza las reflexiones de las ondas y se realiza una comparativa entreambos ejemplos.

5.2 validación experimental de la solución fundamental estratifica-da

En este apartado se pretende validar la implementación de la función de Green delsemiespacio estratificado formulado por Kausel [24, 25] y presentada en el capítuloanterior. Para ello se cuenta con los resultados de una campaña experimental realiza-da para el estudio vibratorio del tramo Palencia-León de la futura LAV de Madrid-Asturias. Aquí se pretende modelizar dicho suelo suponiendo un semiespacio estra-tificado cuyas propiedades dinámicas son las resultantes del estudio, y compararselos resultados teóricos con los experimentales fruto de la caracterización del suelo.Concretamente se comparan las curvas de dispersión calculadas teóricamente con lasobtenidas experimentalmente en el ensayo SASW para validar así la implementacióndel modelo.

101

102 ejemplos numéricas

En el apartado siguiente se explica brevemente cómo se identifican las propiedadesdinámicas del suelo basadas en el ensayo SASW y Sísmica de refracción. Y en el últimoapartado se valida la implementación del modelo teórico, basado en el método híbrido,mediante comparación de las curvas de dispersión experimental.

5.2.1 Identificación de las propiedades dinámicas del suelo

Para identificar las propiedades dinámicas del suelo se suele recurrir a dos méto-dos: el método sísmico de refracción y el análisis sísmico de las ondas superficiales(SASW). Estos métodos requieren de unos ensayos que miden la respuesta del sueloen diferentes puntos tras excitarse por un impulso. Mediante el método sísmico derefracción se puede determinar la velocidad de las ondas P en un suelo estratificado.Éste se basa en medir los tiempos de llegada de las ondas primarias (P) a cada uno delos puntos de medida. Las primeras ondas en registrarse son las ondas P directas paralos puntos más cercanos a la fuente, mientras que para los más alejados suelen ser lasP refractadas. Esto se debe a que normalmente las capas inferiores suelen poseer unasvelocidades de propagación mayores, luego la mayor distancia que deben recorrer lasondas refractadas se compensa con la mayor velocidad de las mismas. De esta forma,se obtiene una curva lineal a trozos que relaciona el tiempo de llegada de las ondasP a cada uno de los puntos de medida. La estratigrafía de la velocidad de las ondasP se calcula al resolver un problema inverso teórico [44]. Por otra parte, mediante elmétodo SASW se puede determinar la velocidad de las ondas S y el amortiguamientode un suelo estratificado. Éste método se basa en el carácter dispersivo de las ondas su-perficiales en un semiespacio estratificado. La respuesta en cada uno de los puntos demedida se transforma al dominio número de onda-frecuencia y se realiza un análisisespectral. Como resultado se obtienen unas curvas de dispersión en las que se repre-senta la velocidad de fase, relacionada con la velocidad de las ondas S, en función dela frecuencia. Nuevamente estas curvas pueden emplearse para obtener la estratigrafíade las ondas S resolviendo un problema inverso [32].En el estudio vibratorio se realizó una caracterización dinámica del suelo mediante

ambos métodos. Para ello se realizó un ensayo in situ que consistió en medir las ace-leraciones en 28 puntos de medida alineados y distanciados cada 3 m. La excitaciónse realizó mediante un martillo calibrado sobre una placa de cimentación cuadradade 300 mm de lado, y las aceleraciones verticales se registraron mediante aceleróme-tros fíjamente unidos a estacas clavadas en el suelo. Se realizaron 100 impactos cuyosresultados fueron promediados en aras a paliar los efectos indeseados del ruido. Delos ensayos se obtuvieron las curvas de dispersión y el tiempo de llega de las curvasde llegada de las ondas P a cada punto. Las propiedades dinámicas fueron halladasresolviendo un problema inverso, y se muestran en la tabla 5.1. Para más detalle sobrelos ensayos, equipos de medida y caracterización dinámica se refiere a Romero [42].

5.2 validación experimental de la solución fundamental estratificada 103

Estrato Espesor [m] cp [m/s] cs [m/s] ξ

1 1.2 1412 189 0.060

2 2.9 1872 189 0.060

3 25.0 1872 495 0.050

4 23.6 1872 440 0.040

5 ∞ 1872 342 0.037

Tabla 5.1: Perfil estratigráfico en el PK 93.850 de la LC Palencia-León [42].

5.2.2 Validación numérica

En este apartado se pretende reproducir teóricamente el ensayo realizado, y obtenerlas curvas de dispersión teóricas del método SASW para compararlas con las obtenidasexperimentalmente. El modelo teórico se realiza en base a la función de Green parael semiespacio estratificado formulado por Kausel y mostrada en el apartado 4. Sesupone un suelo cuyos perfil estratigráfico es el obtenido del estudio. En este puntose debe notar que el coste computacional del método híbrido suele ser muy elevado ydepende de la distancia al punto de medida, o lo que es lo mismo del tiempo máximode evaluación, de las propiedades de los estratos y sobre todo de la profundidad de losmismo. Como puede verse en la tabla 5.1 el tercer y cuarto estrato tienen un elevadoespesor, y sin embargo las propiedades dinámicas son muy similares. En el modelose ha realizado una simplificación, y solo se ha considera 3 estratos. De hecho se harealizado dos modelos: uno en el que se ha considerado los dos primeros estratos y eltercero se ha extendido hasta el infinito, y otro en el que se ha considerado igualmentelos dos primeros estratos bajo el cual se extiende un semiespacio homogéneo con laspropiedades del último estrato de la tabla.

Estrato Espesor [m] cp [m/s] cs [m/s] ξ

Modelo 1

1 1.2 1412 189 0.050

2 2.9 1872 189 0.050

3 ∞ 1872 495 0.050

Modelo 2

1 1.2 1412 189 0.050

2 2.9 1872 189 0.050

3 ∞ 1872 342 0.050

Tabla 5.2: Perfil estratigráfico del modelo teórico

Igualmente que en el ensayo, se ha calculado los desplazamientos verticales para28 puntos de medida separados 3 m, hasta una distancia máxima de 84 de la fuente.Se ha considerado un impulso teórico sintético tipo senoidal de fuerza unitaria quetiene una duración de tp = 0.01 seg y que actúa sobre una área circular de 300 mm dediámetro. Seguidamente la respuesta en aceleración se ha trasformado al dominio de


la frecuencia por aplicación de la trasformada de Fourier y se ha calculado la funciónde transferencia H(r,ω):

H(r,ω) =A(r,ω)

F(r,ω)(5.1)

donde A(r,ω) y F((r,ω) es la trasformada de Fourier a la frecuencia de las acelera-ciones y de la fuerza. En el caso experimental se utiliza la función de transferenciapromediada para los 100 impactos en aras a paliar el ruido de la señal para distanciaslargas entre la fuente y el punto de medida.Para obtener el espectro en el dominio número de onda-frecuencia, a la función

de transferencia se le ha aplicado una trasformada de Hankel que tiene en cuenta lapropia simetría radial del problema tal como se comentó en el apartado 4.4.2:

H(kr ,ω) =∫ ∞

0H(r,ω)J0(kr , r)rdr (5.2)

donde J0(krr) es la función de Bessel de orden cero. Evidentemente la integral debeser truncada para la máxima distancia considerada, lo que tiene cierta influencia enlos resultados si no se escoge una distancia suficientemente grande.Las ondas de Rayleigh está relacionado con los valores máximos del espectro en

el dominio de número de onda-frecuencia [44]. Y las ondas S está relacionada conéstas según un factor que depende del coeficiente ν del suelo. En la Figura 5.1 semuestra el espectro en función de la velocidad de fase kr/ω y de la frecuencia ω delensayo experimental y de los dos modelos teóricos. Puede verse que para frecuenciasinferiores a 15 Hz la resolución del espectro es deficiente para todos los casos. Ademáspuede apreciarse que la tendencia de los máximos del espectro para frecuencias altases similar para los tres casos mostrados.Finalmente las curvas de dispersión se calculan buscando el valor de la velocidad de

fase que maximiza el argumento del espectro para cada frecuencia. Matemáticamentese expresa como:

kr(ω) = Max[H(kr ,ω), kr] (5.3)

Tal que las curvas de dispersión de las ondas superficiales en función de la frecuenciase define como:

cR(ω) =ω

kr(ω)(5.4)

En la Figura 5.2 se muestran una gráfica en la que se superponen las curvas dedispersión del ensayo experimental (línea negra) y de los dos modelos teóricos (lineas

5.2 validación experimental de la solución fundamental estratificada 105

(a) Ensayo experimental (b) Modelo teórico 1

(c) Modelo teórico 2

Figura 5.1: Espectro de onda en el dominio número de onda-frecuencia para el ensayo experimental ypara los modelos teóricos.

gris claro corresponde al modelo 2 y gris oscuro al modelo). Se ha eliminado la bandainferior a 15Hz donde la precisión del espectro es deficiente. Puede apreciarse queefectivamente la tendencia de las curvas de dispersión a frecuencias altas coincidenen los tres casos para y corresponde con la velocidad de fase del estrato superior. Sinembargo, para frecuencias bajas se observa un desfase de los modelos respecto de lacurva experimental, en el que el desfase es mas acentuado para el modelo 1. Se debepresumiblemente a que la velocidad de fase a baja frecuencia depende de las propie-dades de los estratos inferiores. En el primer modelo se realizó una simplificación dela estratigrafía y se homogeneizó los tres últimos tomando las propiedades relativas altercer estrato cuya velocidad de las ondas S es superior a las de los restantes inferiores.


Sin embargo para el modelo 2 se realizó la misma simplificación pero se tomaron laspropiedades del último estrato cuya velocidad. Efectivamente puede verse que paraeste modelo el desfase es menor.

0 20 40 60 80 100100

150

200

250

300

350

400

Frequency [Hz]

Pha

se v

eloc

ity [m

/s]

ExperimentalModelo2Modelo 1

Figura 5.2: Comparación de las curvas de dispersión experimental (línea negra), con las curvas teóricasdel modelo 1 (linea gris claro) y del modelo 2 (línea gris oscura)

Por otro lado, es conveniente estudiar la influencia de la distancia máxima del últimopunto de media en las curvas de dispersión. En la Figura 5.3 se muestra la curva dedispersión experimental (línea negra) frente a las curvas teóricas del modelo 2 parauna separación máxima a la fuente de xmax = 84 m, 66 m, 48 m. Se puede observarque a medida que se disminuye la distancia máxima se acentúa el desfase con laexperimental a baja frecuencia. Se debe a que efectivamente los puntos de medida másalejados son los que tienen en cuenta las propiedades de los estratos inferiores, por loque al perderse esa información afecta a las curvas en la región de baja frecuencia.

En conclusión se tiene que los modelos se aproximan razonablemente bien a lascurvas experimentales . Además en las curvas de dispersión, la velocidad de fase aaltas frecuencias se corresponde con los estratos superiores mientras que a baja con losestratos inferiores. Por tanto la curva teórico es tanto más exacta cuanto evidentementemejor sea el modelo de la estratigrafía y mayor sea la distancia del último punto demedida.

5.3 influencia de la profundidad de la fuente o el receptor

La solución fundamental del semiespacio homogéneo [19], tal como se ha mostrado,se presenta en forma integral, y contiene numerosas problemas numéricos que en lapráctica conllevan una evaluación costosa, sobre todo si debe computarse repetidasveces como sucede en el MEC para estudiar problemas de propagación de ondas en elsuelo. En este caso está muy extendida el uso de la solución fundamental del espaciocompleto debido a la sencillez de la propia formulación [16] y a la relativa facilidad

5.3 influencia de la profundidad de la fuente o el receptor 107

0 20 40 60 80 100100

150

200

250

300

350

400

Frequency [Hz]

Pha

se v

eloc

ity [m

/s]

Figura 5.3: Comparación de las curvas de dispersión experimental (línea negra), con las curvas teóricasdel modelo 2 para las distancias de medida máxima de xmax = 84m (linea gris oscuro), xmax = 66m (lineagris claro) y xmax = 48m (linea gris discontinua)

a la hora de implementarse dentro del MEC. Sin embargo, para tener en cuenta lainfluencia de la superficie libre es necesario mallarla e imponer tensiones nulas sobrela misma, que evidentemente introducen un elevado número de grados de libertad enel modelo y penaliza en la eficiencia del método. Mientras que si se utiliza la soluciónfundamental del semiespacio se satisface implícitamente la condición de superficie li-bre. En consecuencia, parece llegarse a un dilema en cuanto a la implementación deuna u otra solución fundamental en modelos basados en el MEC para problemas deprogación de ondas en el suelo. En este apartado, precisamente se estudia la influenciaque supone la superficie libre del semiespacio, las nuevas ondas que aparecen frutode las reflexiones y su relativa importancia. Concretamente se estudia para qué pro-fundidad tanto de la fuente como del receptor se debe alcanzar para que las ondasreflejadas no tengan apenas repercusión sobre la respuesta.

Para este estudio se ha tomado un suelo cuyas propiedades dinámicas son ρ =1950 kg/m3, cp = 600 m/s, cs = 250 m/s y en el que no se ha tenido en cuenta elamortiguamiento. Se ha considerado que tanto la fuente como el receptor se encuen-tra a la misma cota d variable y separados una distancia horizontal r también variable.Se ha supuesto una excitación vertical tipo escalón unitario y se han registrado losdesplazamientos en la misma dirección. Notar que se podría haber supuesto en sulugar un fuente tipo delta de Dirac, pero a juicio del autor se considera que es menosidónea para este estudio. Ello se debe a que mientras que la respuesta a un impulso sedesvanece con el tiempo, y por tanto los desplazamientos tienden a un permanente nu-lo, para un escalón a parte de poderse estudiar la respuesta dinámica, puede tambiénestudiarse la respuesta estática que es diferente de cero.


La solución a este problema se ha calculado mediante la formulación de Johnsondel semiespacio presentada en el capítulo 3 y mediante las expresiones analíticas delespacio completo que pueden encontrarse en la referencia [16]. Esta última se presentapara un impulso, tanto que aquí ha sido integrada de forma analítica y su expresiónpuede encontrarse en el Anexo A. Ambas soluciones se han comparado a diferentescotas de profundidad y distancias radiales para estudiar la influencia del semiespacio.En la Figura 5.4 se muestra la evolución temporal de los desplazamientos obteni-

das mediantes ambas formulaciones para el caso en el que el receptor se distanciade la fuente 5m horizontalmente. Se muestran nueve curvas relativas a las profun-didades: d = 0, 1, 2, 4, 8, 16, 32, 64, 128 m. Y para cada cota se superponen losdesplazamientos calculados mediante la fórmula de Johnson en línea continua negra,y en discontinua gris oscura los calculados utilizando la fórmula del espacio comple-to. Los desplazamientos y el tiempo se han representado en forma adimensional por(⊓zz = uzzµπr/Z) y por τ = tcs/r, siendo Z la magnitud de la fuerza. Los desplaza-mientos se han evaluado hasta un tiempo igual a dos veces el tiempo de llega de lasondas reflejadas S para el caso más alejado de la superficie. Debido al elevado tiem-po de computo, y en aras a mostrar más claramente los resultados, el eje de abscisasse ha dispuesto en escala logarítmica. Además en línea discontinua gris clara, se hanresaltado los tiempos de llegada de cada tipo de ondas: ondas P directas, S directas,S reflejadas (SS), P reflejadas (PP), P reflejas que cambian de naturaleza (PS) y vice-versa (SP), en las que estas dos últimas son coincidentes. En términos generales, lascurvas presentan dos zonas bien diferenciadas, una región transitoria donde los des-plazamientos varían con el tiempo por la influencia de la llegada de las ondas, y unaregión permanente donde los desplazamientos se estabilizan. Puede observarse, que amedida que la profundidad es mayor, la región transitoria se extiende temporalmentedebido a que las ondas reflejadas tardan más tiempo en alcanzar el punto de evalua-ción, y éstas son cada vez más débiles por la condición de radiación. Sin embargo, lasondas directas siempre llegan al mismo tiempo, llegando las ondas S para τ = 1. Porotro lado, se observa que el desplazamiento transitorio de ambas soluciones son muydiferentes debido a que en la solución del semiespacio se contempla todas las ondas(directas, reflejadas y de superficie llamadas también ondas de Rayleigh), mientras queen la del espacio completo solo existen las ondas directas (P y S). Efectivamente a me-dida que aumenta la profundidad ambas soluciones se van asemejando, siendo cadavez menos influyente las ondas reflejadas. De hecho para una profundidad entorno aa 16 m la influencia de éstas últimas son despreciables, no obstante el desplazamientopermanente o estático de ambas difieren (esta diferencia en el regimen permanenteno se hubieran captado si se hubiera utilizado un impulso). Debe alcanzarse una pro-fundidad de unos 128 m para que los desplazamientos estáticos sean coincidentes yambas soluciones sean por tanto completamente similares.En las Figuras 5.5 y 5.6 se muestran la misma gráfica que la obtenida anteriormen-

te pero para una distancia horizontal de x = 50 m y de x = 100 m respectivamente.


10−1

100

101

0

1

2

4

8

16

32

64

128

Tiempo [−]

Pro

fundid

ad

[m]

Semiespacio

Espacio Completo

P

S

PS

PP

SS

Despla

zam

iento

[--]

0

0.5

Figura 5.4: Comparativa entre la solución del semiespacio y la solución del espacio completo. La fuente yel receptor se encuentran a la misma cota y a 5 m de distancia. Se representa nueve curvas relativas a lasprofundidades: d = 0 m, d = 1 m, d = 2 m, d = 4 m, d = 8 m, d = 16 m, d = 32 m, d = 64 m, d = 128 m.En linea negra continua se representa la evolución temporal de los desplazamientos calculados con laexpresión de Johnson, en linea discontinua gris oscuro la correspondiente al espacio completo, y en grisclaro discontinuo con marcas se representa el tiempo de llegada de las ondas.

Se pretende ahora estudiar cómo influye en ambas soluciones la distancia horizontal.Efectivamente, las observaciones comentadas para la gráfica anterior se vuelven a co-rroborar en éstas, sin embargo observando estas tres figuras comparativamente puedeverse que a medida que se distancian (la fuente respecto del receptor) es mayor lainfluencia de la superficie. Esto es, la influencia de las ondas de Rayleigh y las ondasreflejadas en la región transitoria perduran hasta una profundidad cada vez mayor.E igualmente la diferencia en el desplazamiento estático de ambas soluciones es más


elevada para la misma cota. Con lo cual a medida que se distancia se requiere mayorprofundidad para que las soluciones tengan el mismo grado de acuerdo.

10−1

100

0

1

2

4

8

16

32

64

128

Tiempo [−]

Pro

fundid

ad

[m]

Semiespacio

Espacio Completo

P

S

PS

PP

SS

Despla

zam

iento

[--]

0

0.5


En conclusión, de estas figuras se ha observado que a medida que la profundidades mayor, la influencia de la superficie en la respuesta es menor, y por tanto tiendena converger los desplazamientos calculados mediante las expresiones del semiespaciocon los del espacio completo. Además, a media que se aumenta la distancia horizontalla convergencia se alcanza para un profundidades mayores. De ello se puede deducirque realmente el ratio entre la distancia horizontal y la profundidad es el parámetro


10−1

100

0

1

2

4

8

16

32

64

128

Tiempo [−]

Pro

fundid

ad

[m]

Semiespacio

Espacio Completo

P

S

PS

PP

SS

Despla

zam

iento

[--]

0

0.5


influyente y que determinada el grado de acuerdo. Concretamente para este suelose ha observado que para un ratio superior a r/d = 5/16 la influencia de las ondassuperficiales y reflejadas en la respuesta transitoria es despreciable, y para un ratiosuperior a r/d = 5/128 los desplazamientos estáticos son prácticamente coincidentes.


5.4 influencia del amortiguamiento del suelo

La formulación presentada por Johnson para evaluar los desplazamientos del semi-espacio es válida para un semiespacio elástico en asuencia de atenuación. Sin embargorealmente el suelo suele poseer un comportamiento inelástico disipativo asociado a laatenuación de los desplazamientos con el tiempo. En este sentido, las expresiones pro-puestas por Kausel[24, 25], y presentadas en el apartado 4, permiten la posibilidad deintroducirle el efecto de un amortiguamiento a la respuesta del semiespacio. En esteapartado se muestra cómo se atenúan los desplazamientos al introducirle amortigua-miento. Esta atenuación no debe ser confundida con la que implícitamente requiere lacondición de radiación de Sommerfeld [16], y que asegura que las ondas se propagandesde el emisor hasta el infinito, tendiendo los desplazamientos asintóticamente a cerocon la distancia.Para ello, se ha supuesto un suelo con las propiedades dinámicas: ρ = 1950 kg/m3, cp =

467.70 m/s, cs = 250 m/s, excitado mediante un impulso tipo senoidal sobre la superfi-cie en la dirección vertical. Se han calculado los desplazamientos en la misma direcciónmediante las expresiones presentadas por Kausel [24, 25] del semiespacio (y mostradasen el apartado 4) en un punto superficial para diferentes casos de amortiguamiento.En la Figura 5.7 se muestra la evolución temporal de los desplazamientos para un

punto situado a 6 m de distancia de la fuente. Se dibujan once curvas relativas adiferentes coeficientes de amortiguamiento comprendidos entre ξ = [0, 0.1]. Se ha eva-luado hasta un tiempo igual a dos veces el tiempo de llegada de las ondas S. Losdesplazamientos y el tiempo se han adimensionalizado por Uzz = uzzµπx/Z y porτ = tcs/x, siendo Z la magnitud de la fuerza y x la separación del receptor respectoa la fuente. La primera observación clara es que los desplazamientos se atenúan a me-dida que se aumenta el amortiguamiento. Además notar que las curvas se suavizan,es decir, la variación de los desplazamientos respecto del tiempo es menos acusada alaumentar el amortiguamiento, o lo que es lo mismo la magnitud de la velocidad delos desplazamientos se atenúa. En este sentido, fíjese que la respuesta de un impul-so corresponde con la derivada temporal de los desplazamientos obtenidos para unaexcitación tipo escalón, esto es, corresponde con la velocidad de los desplazamientospara esta última excitación. En la Figura 5.8 se muestra la respuesta a una excitaciónescalón obtenida mediante integración de la respuesta al impulso. Observar que du-rante el régimen transitorio, donde los desplazamientos varía con el tiempo por elefecto de llegada de cada una de las ondas, se produce una atenuación de los mismoy se suavizan las curvas. Sin embargo el desplazamiento permanente permanece inal-terable y alcanza siempre el mismo valor adimensional. Esto quiere decir que el efectode introducir un amortiguamiento hace que se atenúe la magnitud de la velocidad dedesplazamiento (se opone al cambio de los desplazamientos), lo que conlleva a que losdesplazamientos en el regimen transitorio sea menor, pero el desplazamiento estáticosea siempre el mismo.

5.4 influencia del amortiguamiento del suelo 113

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0

0.05

0.1

Tiempo [−]

Fac

tor

de a

mor

tigua

mie

nto

criti

co ξ

Figura 5.7: Evolución de los desplazamientos verticales debido a una excitación impulsiva vertical paradiferentes factores de amortiguamiento. La fuente y el receptor se encuentra sobre la superficie y a 6 mde distancia. Se representa once curvas relativas a los factores de amortiguamiento: ξ = 0, ξ = 0.01, ξ =0.02, ξ = 0.03, ξ = 0.04, ξ = 0.05, ξ = 0.06, ξ = 0.07, ξ = 0.08, ξ = 0.09, ξ = 0.1.

Por otro lado, debe comentarse que la atenuación de la velocidad por el efecto delamortiguamiento material no tiene nada que ver con la atenuación a la condición deradiación de Sommerfeld que debe satisfacer implícitamente las ecuaciones de propa-gación de ondas en un semiespacio infinito. Esta condición asegura que las ondas sepropagen desde el emisor hacia el infinito. A medida que el frente de onda se aleja de


0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0

0.05

0.1

Tiempo [−]

Fac

tor

de a

mor

tigua

mie

nto

criti

co ξ

Figura 5.8: Evolución de los desplazamientos verticales debido a una excitación tipo escalón vertical paradiferentes factores de amortiguamiento. La fuente y el receptor se encuentra sobre la superficie y a 6 mde distancia. Se representa once curvas relativas a los factores de amortiguamiento: ξ = 0, ξ = 0.01, ξ =0.02, ξ = 0.03, ξ = 0.04, ξ = 0.05, ξ = 0.06, ξ = 0.07, ξ = 0.08, ξ = 0.09, ξ = 0.1.

la fuente la energía de las ondas se reparte entre una superficie mayor, por lo que dis-minuye la velocidad de desplazamiento de de las partículas. Tal como se expresó en elapartado 2, la condición de radiación implica que la velocidad de los desplazamientosdecaigan inversamente proporcional a la distancia a la fuente. En otras palabras, elamortiguamiento material está asociado a una disipación de energía mientras que la

5.5 suelo compuesto por un estrato blando sobre semiespacio rocoso 115

condición de radiación no. En la Figura 5.9 se representa la evolución de los desplaza-mientos superficiales con el tiempo (eje x) y con la separación a la fuente (eje y) debidoa una excitación superficial vertical impulsiva sin amortiguamiento. Los desplazamien-tos se ha adimensionalizado por Uzz = uzzµπxmax/Z, el tiempo por τ = tcs/xmax, yla distancia a la fuente por X = x/xmax, siendo xmax la distancia máxima evaluada.De nuevo resaltar que los desplazamientos debido al impulso corresponden con lasvelocidades producidas por una excitación escalón. Efectivamente puede verse cómola velocidad decae inversamente proporcional a la distancia debido a la carga escalón.

Figura 5.9: Evolución de los desplazamientos verticales con respecto al tiempo (eje x) y en función de ladistancia a la fuente (eje y) debido a una carga impulsiva. La fuente y el receptor se encuentra sobre lasuperficie. Se ha considerado un amortiguamiento material nulo.

5.5 suelo compuesto por un estrato blando sobre semiespacio rocoso

En este apartado se estudia la propagación de ondas en un suelo compuesto por unestrato blando que apoya sobre un semiespacio rocoso. Se calcula la evolución de losdesplazamientos superficiales con el tiempo mediante el método híbrido propuesto


por Kausel y presentado en en el apartado 4, y se analiza las diferentes ondas reflecta-das que surgen.Concretamente se supone un estrato de 5 m de espesor con las propiedades: ρ =

1500 kg/m3, cp = 336.75 m/s, cs = 180 m/s, ξ = 0.005 , sobe un estrato rocoso semi-infinito cuyas propiedades son: ρ = 1950 kg/m3, cp = 1216 m/s, cs = 650 m/sξ =

0.005. Se supone que el semiespacio es excitado mediante una carga vertical impul-siva tipo senoidal sobre la superficie (que puede asemejarse con un martillo de im-pacto), y se registran los desplazamientos verticales superficiales para las distanciasde x = 2, 4, 6, 8, 10 m. En la Figura 5.10 se muestra cinco curvas que represen-tan la evolución temporal de los desplazamientos para las diferentes distancias con-sideradas. Los desplazamientos y el tiempo se muestran de forma adimensional porUzz = uzzµπxmax/Z y por τ = tcs/xmax. En línea discontinua gris se muestra el tiempode llegada de cada una de las ondas que aparecen, excepto las ondas de Rayleigh queevidentemente coincide con los desplazamientos máximos y solo se propagan en lasuperficie.En primer lugar se debe observar la complejidad de la respuesta donde aparecen

muchos tipos de ondas. Inicialmente, al aplicar la excitación se producen tres tiposde ondas, las primarias (P), las secundarias o transversales (S) y las de superficie (On-das de Rayleigh). Éstas últimas solo se propagan por la superficie y no pueden sufrirreflexiones. Sin embargo las ondas P y S se propagan dentro del medio (a diferentevelocidad evidentemente), de forma que al encontrarse con otro medio de propiedadesdiferentes se produce el fenómeno de reflexión, en el que las ondas rebotan y regresana la superficie, y el fenómeno de refracción en el que parte de las ondas se propaganpor el otro medio. Evidentemente, las ondas a medida que se distancian de la fuen-te se van atenuado por el efecto del amortiguamiento material y por la condición deradiación, por ello las reflexiones de las ondas son cada vez menos importantes talcomo se observa. En este caso particular, para las distancias que se han considerado(x = [2 m, 10 m]) y para la ventana temporal donde se han calculado los desplaza-miento (2 veces el tiempo de llega de las ondas S para el caso de la máxima distancia),solo pueden registrarse como máximo tres reflexiones, es decir, las ondas emitida porla fuente pueden rebotar en el estrato rocoso, volver y rebotar en la superficie, y denuevo rebotar en el estrato para volver a la superficie. Además, como se explicó en elapartado 3, las reflexiones pueden ser de dos tipos, uno en la que la naturaleza de lasondas no cambia ,y otro en la que cambia P ↔ S. En el primer tipo el ángulo de inci-dencia coincide con el de reflexión, mientras que en el otro los ángulos se relacionansegún la Ecuación (3.97). Concretamente para este ejemplo existen 23 tipos de ondasde las cuales 20 de ellas son ondas reflejadas. En la tabla 5.3 se muestran todos lostipos de ondas:donde cada letra representa la naturaleza de la onda que se propaga (P o S) hasta

que se produce la reflexión y que puede cambiar la naturaleza de la misma , es decir,las ondas PS corresponde con las ondas que originariamente son primarias pero al

5.5 suelo compuesto por un estrato blando sobre semiespacio rocoso 117

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

2

4

6

8

10

Tiempo [−]

Dis

tanc

ia [m

]

SemiespacioPSPSPPSSPSPSPSSSSPPPPPPP

Figura 5.10: Evolución temporal de los desplazamientos verticales de un semiespacio estratifica com-puesto por un estrato blando de 5 m de espesor sobre un suelo rocoso debido a una excitación impulsivatipo senoidal vertical. La fuente y el receptor se encuentra sobre la superficie y separados una distanciavariable (x = 2, 4, 6, 8, 10 m).

Número de reflexiones tipo de ondas

0 P S Rayleigh

1 PP SS PS SP

3PPPP PPSS PPSP PPPS PSPP PSSS PSSP PSPS

SPPP SPSS SPSP SPPS SSPP SSSS SSSP SSPS

Tabla 5.3: Tipo de ondas y reflexiones

reflejarse cambian a ondas secundarias. No obstante, para este caso, donde la fuentey el receptor se encuentran en la superficie, el tiempo de llegada de ciertas ondas


coinciden, como puede verse en la figura. En la tabla 5.4 se muestra la coincidencia enel tiempo de llegada entre las distintas ondas reflejadas.

Tiempo representativo tipo de ondas equivalentes

tPS PS SP

tPSPS PSPS SPSP PPSS SSPP PSSP SPPS

tPSSS PSSS SPSS SSPS SSSP

tSPPP SPPP PSPP PPSP PPPS

Tabla 5.4: Equivalencias entre los tiempos de llegadas de las ondas reflejadas

Por otro lado, en la Figura 5.11 se muestra una comparativa de la evolución temporalde los desplazamientos calculados para este suelo respecto a los calculados suponien-do que el estrato blando es libre, es decir, no se apoya sobre el suelo rocoso, y parauna distancia de 4 m de la fuente. Los desplazamiento y tiempo se muestra de formaadimensional. Se puede observar en la figura que las ondas directas producen el mis-mo campo de desplazamiento en ambos modelos. Efectivamente estas ondas no tienenen cuenta la influencia de la superficie inferior, siempre y cuando ésta se localice su-ficientemente profunda. En cuanto a las ondas reflejadas obviamente llegan al mismotiempo pero el valor de los mismos difieren. Es más, si se observa detenidamente sepercata que la magnitud de los desplazamientos en valor absoluto son similares, peropara el caso del estrato libre se invierte la dirección de los desplazamientos tras lareflexión. En linea gris claro de punto y raya se ha invertido la dirección de los des-plazamientos de las ondas reflejadas de la placa, y puede verse que efectivamente sonmuy simulares respecto a las del suelo en cuestión. Ello indica que las ondas se propa-gan preferentemente por el estrato blando gracias al efecto de la reflexión, siendo portanto bajo la proporción de ondas refractadas al medio rocoso.

5.6 suelo con inclusiones rocosas

Por último en este apartado se estudia la propagación de ondas en un suelo blandoen el que se le intercala un fino estrato rocoso. Las propiedades de los estrados sonlas mismas que las tomadas para el apartado anterior. El estrato rocoso se localiza a5 m de la superficie y cuenta con un espesor de 0.5 m. Igulamente se ha consideradouna excitación impulsiva tipo senoidal sobre la superficie en dirección vertical, y sehan calculado los desplazamientos superficiales en la misma dirección mediante elmétodo híbrido formulado por Kausel.En la Figura 5.12 se representa nuevamente la evolución de los desplazamientos

con el tiempo para diferentes distancias del receptor respecto a la fuente. Los des-plazamientos se muestran adimensional frente al tiempo también adimensional. Enlinea discontinua gris claro se marca el tiempo de llegada de las ondas directas y lasreflejadas en la primera interfase. En rasgos generales puede verse que el campo de

5.6 suelo con inclusiones rocosas 119

0 0.5 1 1.5 2−500

0

500

1000

Tiempo [−]

Des

plaz

amie

nto

[m]

Suelo estratificadoEstrato libre−libreInversa de las reflexiones

Figura 5.11: Comparativa de la evolución temporal de los desplazamientos de un semiespacio estratifi-cado (linea continua) respecto a los desplazamiento de un estrato libre-libre (linea discontinua gris). Elsemiespacio está compuesto por un estrato blando de 5 m de espesor sobre un suelo rocoso, mientrasque la placa libre tiene el mismo espesor y propiedades que el estrato blando. La excitación es superficialvertical impulsiva de tipo seniodal y se registra los desplazamientos verticales superficiales a 4 m dedistancia. En linea de punto y raya gris claro se representa la inversa de los desplazamientos debido alas reflexiones en la placa.

desplazamiento resultante es muy simular al calculado anteriormente donde se ha su-puesto que el estrato blando se apoya sobre un semiespacio rocoso. Efectivamente laperturbación ocasionada por las ondas directas (P, S y Superficiales) son idénticas, eiguales a su vez a las correspondiente a un semiespacio homogéneo. Sin embargo apriori cabe esperarse que las ondas reflejadas sean diferentes. De hecho en este ejem-plo existen dos superficies de reflexión relativa a las interfases del estrato rocoso, conlo cual pueden existir una mayor variedad de ondas reflejadas. En la Figura 5.13 semuestra una comparativa con la evolución de los desplazamientos del ejemplo ante-rior para una separación de 4 m. Curiosamente no parecen observarse nuevos tiposde ondas reflejadas, y además la magnitud de los desplazamientos provocados porlas ondas reflejadas son muy simulares. Ello parece indicar que la mayor parte de laenergía de las ondas se refleja en el estrato rocoso, suponiendo una fuerte barrera, yuna pequeña parte imperceptible, se transmite al estrato rocoso.


0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

2

4

6

8

10

Tiempo [−]

Dis

tanc

ia [m

]

SemiespacioPSPSPPSSPSPSPSSSSPPPPPPPP

Figura 5.12: Evolución temporal de los desplazamientos verticales de un semiespacio estratifica debidoa una excitación impulsiva tipo senoidal vertical. El semiespacio está compuesto suelo blando en el quese le intercala un delgado estrato rocoso de 0.5 m de espesor y a 5 m de profundidad. La fuente y elreceptor se encuentra sobre la superficie y separados una distancia variable (x = 2, 4, 6, 8, 10 m).

5.6 suelo con inclusiones rocosas 121

P S PP PS SS,PPPP SPPP PSPS SSSP−200

−100

0

100

200

300

400

500

Tiempo [−]

Des

plaz

amie

nto

[m]

Suelo blando con inclusion de estrato rocosoSuelo blando sobre roca

Figura 5.13: Comparativa de la evolución temporal de los desplazamientos de un semiespacio estratifi-cado compuesto por un suelo blando en el que se le intercala un estrato rocoso de 0.5 m de espesor y a5 m de profundidad (linea continua), respecto de un suelo blando de 0.5 m de espesor que apoya sobresemiespacio rocoso infinito. Las propiedades de los estratos son idénticas. La excitación es superficialvertical impulsiva de tipo seniodal y se registra los desplazamientos verticales superficiales a 4 m dedistancia. En linea de punto y raya gris claro se representa la inversa de los desplazamientos debido alas reflexiones en la placa.

6CONCLUS IONES

6.1 conclusiones

En este trabajo se ha estudiado el problema de Lamb [27] sobre propagación deondas en un semiespacio excitado por una carga dinámica. Para ello, el texto se haorganizado en 3 bloques. En el primero se ha presentado la solución fundamental delsemiespacio homogéneo de Johnson. La solución que propone se expresa en formaintegral y está formada genéricamente por seis términos relativos a las ondas direc-tas P y S, y a las ondas reflejadas en la superficie PP, SS PS y SP. Los dos primerostérminos, desde el punto de vista de la implementación y evaluación numérica, nopresenta problema alguno, de hecho Johnson resuelve las integrales analíticamente.Sin embargo los núcleos integrales del resto de términos presentan singularidades, loque dificulta su evaluación numérica. Se ha mostrado que las singularidades débilespueden ser eliminadas mediante un cambio de variables. Mientras que la singularidaddebido a las ondas de Rayleigh puede ser evitado tomando un camino de integracióncomplejo, lo que supone un ahorro considerable de tiempo de cálculo y mejora dela precisión cuando tanto la fuente como el receptor tienden a la superficie. Por otrolado, se ha presentado la solución para el caso particular en el que bien la fuente obien el receptor se encuentran sobre la superficie. Se ha mostrado que la solución ala que se llega es una simplificación de la solución general y su evaluación es menoscostoso computacionalmente. La solución de Johnson se ha programado en Matlab, seha comparado con la solución explícita de Pekeris, y se ha comprobado que se llega alos mismo resultados numéricos.

En el segundo bloque se ha presentado la formulación del semiespacio estratificadopropuesto por Park y Kausel basada en el método híbrido. Este método se planteaen el dominio número de onda-tiempo y combina dos soluciones, la respuesta delsemiespacio homogéneo y la solución de una placa horizontalmente estratifica. Se hamostrado que la formulación del semiespacio homogéneo no es explícita y requie-re resolver integrales impropias fuertemente oscilatorias. Estas integrales pueden serresueltas dividiéndola en partes y aproximádolas por determinadas funciones asintóti-cas tal como propone Kausel, sin embargo el coste computacional es muy alto. En estesentido se ha mostrado que las integrales pueden ser precomputadas y almacenadasen tablas lo que reduce el tiempo de cálculo. Por otro lado, la solución de la placaestratificada se basa en el TLM que conlleva a una discretización parcial del medio enla dirección estratificada y una superposición modal. Se ha mostrado que medianteeste método se obtiene una solución aproximada tanto más exacta cuanto mejor sea la

123

124 conclusiones

discretización y mayor sea el número de modos tenidos en cuenta. Ambos problemasse acoplan resolviendo un sistema de ecuaciones recursivo para cada paso de tiempo,lo que introduce adicionalmente una discretización temporal. Se ha mostrado que estemétodo no es incondicionalmente estable y depende del paso temporal para resolverel sistema de ecuaciones y del numero de elementos tomados por longitud de onda.Por último, en el tercero bloque se ha validado la solución del semiespacio estra-

tificado con resultados experimentales. Se ha comparado la curva de dispersión ex-perimental del ensayo SASW con la obtenida numéricamente y se ha concluido queambas curvas presentan un grado de coincidencia razonable. Además se ha realiza-do una comparativa de la solución del espacio completo con la del semiespacio y seha mostrado que para cierta relación entre la separación horizontal fuente-receptor yprofundidad, la influencia de la superficie tiende a desaparecer.

B IBL IOGRAF ÍA

[1] Lecons sur la Thórie Mathématique de L’elasticité des Corps Solides. PhD thesis, Paris,1852.

[2] J.D. Achenbach. Wave propagation in elastic solids. 1973.

[3] R. Aspel and E. Luco. On the Green’s function for a layared half-space. Part I.Bulletin of the Seismological Society of America, 73:909–929, 1983.

[4] R. Aspel and E. Luco. On the Green’s function for a layared half-space. Part II.Bulletin of the Seismological Society of America, 73:931–951, 1983.

[5] M. Biot. Theory of propagation of elastic waves in a fluid-satured porous solid. I.Low-frecuency range. Journal of Acoustical Society of America, 28:168–178, 1956.

[6] M. Biot. Theory of propagation of elastic waves in a fluid-satured porous solid. II.High-frecuency range. Journal of Acoustical Society of America, 28:179–191, 1956.

[7] M. Bouchon. A simple method to calculate Green’s function for elastic layeredmedia . Bulletin of the Seismological Socmty of America, 71:959–971, 1981.

[8] L. Cagniard. Reflection and Refraction of Progressive Seismic Waves. Mc-Graw Hill,New York, 1962.

[9] C.C. Chao. Dynamical response of an elastic half-space to tangential surfaceloadings. Journal of Applied Mechanics, Transactions of the ASME, 27(3):559–567,1960.

[10] M. Crisfield. An arch-length method including line searches and accelerations .International Journal for Numerical Methods, pages 1269–1289, 1983.

[11] A.T. de Hoop. A modification of Cagniard’s method for solving seismic pulseproblems. Applied Scientific Research, Section B, 8(1):349–356, 1960.

[12] A.T. de Hoop. Theoretical determination of the surface motion of a uniformelastic hald-space produced by a dilatational, impulsive, point source. Colloques

Internationauax du Centre National de la Recherche Scientifique, pages 21–32, 1961.

[13] M. Doblaré. Fundamentos de la elasticidad lineal . SINTESIS, 1998.

[14] J. Domínguez. Boundary elements in dynamics. Elsevier, 1993.

125

126 bibliografía

[15] J. Domínguez and R. Gallego. Time domain boundary elements: a comparative study,.in Advances in boundary elements vol. 3, Brebbia CA. Connor JJ. (eds), Springer,1989.

[16] A.C. Eringen and E.S. Suhubi. Elastodynamics, Volume 2, Linear theory. AcademicPress, New York, USA, 1975.

[17] Y. Hisada. An efficient method for computing Green’s functions for a layaredhalf-space with sources and recievers at close depths. Bulletin of the Seismological

Society of America, 84:1456–1472, 1994.

[18] Y. Hisada. An efficient method for computing Green’s functions for a layaredhalf-space with sources and recievers at close depths (Part 2). Bulletin of the Seis-

mological Society of America, 84:1456–1472, 1994.

[19] L.R. Johnson. Green’s Function for Lamb’s Problem. Geophysical Journal of the

Royal Astronomical Society, 37(1):99–131, 1974.

[20] E. Kausel. An explicit solution for the Green’s functions for dynamic loads inlayered media. MIT Research Report R81-13.

[21] E. Kausel. Wave propagation in anisotropic layered media. International Journal

for Numerical Methods, 23:1567–1578, 1986.

[22] E. Kausel. Thin Layared Method: Formulation in the time domain. International

Journal for Numerical Methods in Engineering, 37:927–941, 1994.

[23] E. Kausel. Fundamental solutions in elastodynamics: a compendium. Cambridge Uni-versity Press, New York, 2006.

[24] E. Kausel and J. Park. Response of layered half-space obtained directly in the timedomain, part I: SH- Sources. Bulletin of the Seismological Society of America, 96(5):1795–1809, 2006.

[25] E. Kausel and J. Park. Response of layered half-space obtained directly in the timedomain, part II: SV-P and three-dimensional sources. Bulletin of the Seismological

Society of America, 96(5):1810–1826, 2006.

[26] E. Kausel and J.M. Roësset. Stiffness matrices for layered soils. Bulletin of the

Seismological Society of America, 71(6):1743–1761, 1981.

[27] H. Lamb. On the propagation of tremors over the surface of an elastic solid.Philosophical Transactions of the Royal Society of London, Serie A, CCIII(1):1–42, 1904.

[28] J. Lysmer. Lumped mass method for Rayleigh waves. Bulletin of the Seismological

Society of America, 60:89–104, 1970.

bibliografía 127

[29] WJ. Mansur. Time-Stepping technique to solve wave propagation problems using the

Boundary Element Method. PhD thesis, 1983.

[30] M. Marrero and J. Domínguez. Numerical behavior of time domain BEM forthree-dimensional transient elastodynamic problems. Engineering Analysis with

Boundary Elements, 27(1):39–48, 2003.

[31] H.M. Mooney. Some numerical solutions for Lamb’s problem. Bulletin of the

Seismological Society of America, 64(2):473–491, 1974.

[32] Nazarian, Soheil, Yuan, Deren, Baker, and R. Mark. Automation of spectral analy-sis of surface waves method. Number 1213, pages 88–100, 1994.

[33] F. París. Teoría de la elasticidad. Grupo de Elasticidad y Resistencia de Material,1996.

[34] F. Paris and J. Cañas. Boundary element methods. Fundamentals and Applica-tions. 1997.

[35] J. Park. Wave Motion in finite and infinite media using the Tin-Layer Method. PhDthesis, February 2002.

[36] J. Park and E. Kausel. Greens Function in the wavenumber-time domain for loadson the surface of a homogeneous half-sapce. 16th Engineerin Mechanics Conference.

ASCE, July 2003.

[37] J. Park and E. Kausel. Impulse Response of Elastic Half-Space in the WaveNumber-Time Domine. Journal of Engineering Mechanics, 130:1211–1222, 2004.

[38] C.L. Pekeris. The seismic surface pulse. Proceedings of the National Academy of

Sciences of the United States of America, 41:469–480, 1955.

[39] W.L. Pilant. Elastic Waves in the earth. 6, 1979.

[40] S. Poisson. Mémoire sur Léquilibre et le Mouvement des Corps Élastiques. Mé-

moires de lAcadémie des Sciences, 8:357–570, 1829.

[41] S. Poisson. Addiction au Mémoiresur Léquilibre et le Mouvement des CorpsÉlastiques. Mémoires de lAcadémie des Sciences, 8:623, 1829.

[42] A. Romero, P. Galvín, and J. Domínguez. Short span bridges dynamic behaviouraccount for the vehicle-track-structure-soil dynamic interaction [Comportamientodinámico de viaductos cortos considerando la interacción vehículo-vía-estructura-suelo]. Revista Internacional de Metodos Numericos para Calculo y Diseno en Ingenieria,28(1):55–63, 2012.

128 bibliografía

[43] M. Schevenels. The Impact of Uncertain Dynamic Soil Characteristics on the Predi-

tion of Ground Vibration. PhD thesis, FAKLTEIT INGENIEURSWTWNSCHAPPEN,May 2007.

[44] M. Schevenels, G. Lombaert, G. Degrande, and S. François. A probabilistic assess-ment of resolution in the SASW test and its impact on the prediction of groundvibrations. Geophysical Journal International, 172(1):262–275, 2008.

[45] L.F. Shampine. Vectorized adaptive quadrature in MATLAB. Journal of Compu-

tational and Applied Mathematics, 211(2):131–140, 2008.

[46] G. Stokes. On the Dynamical Theory of Diffraction. Transactions of the Cambridge

Philosophical Society, 9:1–62, 1849.

[47] W. Thomson. Transmission of elastic wavs through a stratified soil medum. Jour-nal of Applied Physics, 21:89–93, 1950.

[48] T. Triantafyllidis. 3-D time domain BEM using half-space Green’s functions. En-gineering Analysis with Boundary Elements, 8(3):115–124, 1991.

[49] P. Xu and A.K. Mal. Calculation of the in-plane Green’s functions for a layaredviscoelastic solid. Bulletin of the Seismological Society of America, 77:1823–1837, 1987.

ESTUDIO NUMÉRICO DE LA PROPAGACIÓN DE...

Documents

Transcript of ESTUDIO NUMÉRICO DE LA PROPAGACIÓN DE...