Ondas elásticas y acoplamiento fluido-sólido

Referencias: Virieux 1986 Moczo et al. 2007 Cremer et al. 2005 Brekhovskikh y Godin 1990 van Vossen et al. 2002

Un fluido transporta una onda; un sólido, una familia de ellas. Esta guía cubre el compañero elástico de la simulación de ondas FDTD 2D: elastic_fdtd_simulation integra el sistema velocidad-esfuerzo en deformación plana de Virieux (1986) sobre la misma malla escalonada que el esquema acústico, y con él la onda P de compresión, la onda S de cizalla, la onda de Rayleigh de una superficie libre y la onda de Scholte de un contacto fluido-sólido emergen de la misma actualización. Como un fluido es simplemente el límite sin cizalla de las ecuaciones elásticas, un par de mapas de material acopla agua, aire, sedimento y acero en una sola ejecución: la reflexión oblicua con conversión de modo, las ondas de interfase del lecho marino y la transmisión de placas sumergidas salen de los mapas sin más, cada una anclada a su forma cerrada exacta por la batería de validación.

La guía supone el vocabulario del esquema acústico (fuentes, sondas, máscaras de obstáculos, esponjas y la cota de Courant); lee primero la simulación de ondas FDTD 2D si todo ello te resulta nuevo. Las formas cerradas que este esquema contrasta viven en las guías de edificación, de vibración y de acústica submarina: la ley de masas de la predicción del aislamiento de paneles, las ondas de flexión de la transmisión de onda de flexión en uniones de placas, y el lecho marino con cizalla que el modelo de reflexión de fondo fluido de la propagación submarina del sonido no puede representar.

1. El esquema velocidad-esfuerzo P-SV

Un sólido transporta más de una onda. Junto a la onda P de compresión a $c_{P} = (λ + 2 μ) / ρ$ se propaga una onda S de cizalla a $c_{S} = μ / ρ$ , y toda superficie libre de tracción guía una onda de Rayleigh justo por debajo de $c_{S}$ . elastic_fdtd_simulation integra el sistema velocidad-esfuerzo 2D en deformación plana de Virieux (1986, Ec. 2),

ρ \frac{\partial v}{\partial t} = \nabla \cdot τ, \frac{\partial τ}{\partial t} = λ (\nabla \cdot v) I + μ (\nabla v + \nabla v^{T}),

sobre la misma disposición escalonada que el esquema acústico: los esfuerzos normales txx/tyy comparten los centros de celda donde vive la presión acústica, las velocidades ocupan las caras y el esfuerzo de cizalla txy las esquinas (la celda totalmente escalonada de Virieux, desplazada media celda). La API replica fdtd_simulation pieza a pieza: mapas de velocidades c_p y c_s (más rho) en lugar de c, ExplosionSource (una inyección isótropa de esfuerzo) y ForceSource (una fuerza de cuerpo direccional) como fuentes, sondas que registran p, vx o vy, la misma maquinaria de esponjas y obstáculos, y un ElasticFDTDResult congelado con el mismo .plot(). Lo que la ejecución anota viaja en un ElasticRecording (las celdas de sonda y sus campos, la cadencia de instantáneas y su campo) y cómo terminan los bordes del dominio en un ElasticBoundaries (los lados y el grosor de la esponja).

Dos números gobiernan cada ejecución, y leen extremos opuestos de los mapas de material. Estabilidad: el paso temporal sigue a la onda más rápida del dominio, $c_{P}^{m a x}$ , exactamente igual que en el esquema acústico (Virieux, Ecs. 6-7). Resolución: el tamaño de celda sigue a la más lenta, $Δ x \leq c_{min} / (10 f)$ , donde $c_{min}$ es la velocidad de cizalla no nula más lenta de entre las celdas sólidas, o una velocidad $P$ de fluido si alguna región fluida es todavía más lenta. En un mapa completamente sólido la longitud de onda de cizalla es siempre la más corta; en un mapa mixto no tiene por qué serlo, como con agua a 1480 m/s sobre un sedimento cuya velocidad de cizalla es de 2000 m/s.

Esos dos extremos son lo que cuesta un dominio de contraste alto. Δx lo fija el medio más lento y Δt el más rápido, así que el número de pasos que hacen falta para cubrir un tiempo físico dado crece con la razón entre las dos velocidades: un dominio de agua sobre acero resuelto para la longitud de onda del agua sigue teniendo que avanzar a los 5900 m/s del acero, aproximadamente cuatro veces más pasos que ese mismo dominio lleno solo de agua. Las regiones lentas trabajan entonces a un número de Courant local mucho menor, donde el error de dispersión del esquema está más cerca de su cota en el peor caso, así que su resolución efectiva es peor de lo que sugiere el recuento de celdas: por eso la regla de las diez celdas es un suelo y no un objetivo en los problemas de contraste alto. Cuando solo uno de los medios lleva la física de interés, sustituir el otro por una condición de contorno (un borde de impedancia o una esponja) suele salir mucho más barato que mallar los dos.

import numpy as np
from phonometry import simulation

# Un bloque de aluminio golpeado por una pequeña explosión: el frente P
# alcanza dos sondas a 100 y 220 celdas; c_P sale del retardo diferencial.
w = 8e-6
res = simulation.elastic_fdtd_simulation(
    6320.0, 3130.0, 0.002, 1.1e-4, rho=2700.0, shape=(501, 501),
    sources=[simulation.ExplosionSource(
        ix=250, iy=250,
        waveform=simulation.GaussianPulse(0, 0, width=w).value)],
    recording=simulation.ElasticRecording(
        probes=[(350, 250), (470, 250)], probe_fields=("p",)),
)
t = res.times[1:]
t1, t2 = (t[np.abs(res.signals[k, 0, 1:]).argmax()] for k in range(2))
print(round(0.24 / (t2 - t1)))    # 6316  c_P = sqrt((lambda+2mu)/rho) = 6320

2. Fluidos dentro del esquema elástico

Poner c_s = 0 marca una celda de fluido: el módulo de cizalla se anula, el sistema degenera en las ecuaciones acústicas y, con un mapa uniforme de c_s = 0, el esquema elástico reproduce el acústico bit a bit (la batería de regresión comprueba la igualdad exacta). Esa equivalencia es una garantía de corrección, no una recomendación: este esquema lleva tres campos de esfuerzo donde el acústico lleva una sola presión, así que aquí un problema enteramente fluido cuesta aproximadamente el doble de memoria y de tiempo de avance. Usa el esquema acústico salvo que alguna parte del dominio soporte cizalla.

Una columna de agua sobre un semiespacio de acero son, por tanto, solo dos bandas de los mapas de material: con la densidad promediada aritméticamente a las caras y el módulo de cizalla armónicamente a las esquinas (Moczo et al. 2007, Ecs. 7.37-7.39), la continuidad de tracción de toda interfase interna, contactos fluido-sólido incluidos, emerge de los mapas sin más. La batería de validación mide la reflexión a incidencia normal de una interfase agua-acero con un error inferior al 2 % (típicamente una fracción de punto porcentual) respecto a $(Z_{2} - Z_{1}) / (Z_{2} + Z_{1})$ y recupera la ley de masas a incidencia normal de una placa de acero de 3 mm sumergida con un error inferior a 0,3 dB respecto a mass_law_transmission_loss.

Dónde está realmente el contacto, y cómo comprobarlo. Una región pintada desde la fila $i$ hacia abajo sitúa la interfase física en el plano de caras $y = i Δ x$ , porque es ahí donde la densidad se promedia a las caras y el módulo de cizalla a las esquinas. Las distancias que se citen «desde el contacto» hay que medirlas desde ese plano, y una sonda tiene que quedar varias celdas apartada de él para que nunca muestree las celdas promediadas. Comprueba un mapa pintado antes de ejecutarlo, con plt.imshow(sim.c_s) y los bordes de las regiones dibujados encima, o simplemente verificando el índice de fila de la primera fila sólida, porque un error de una fila desplaza todas las llegadas en $Δ x / c$ y cuesta una ejecución entera. No hay contrapartida elástica de FDTD2D.plot_geometry, así que esta comprobación corre por cuenta del lector. Merece la pena conocer de antemano un fallo más: toda celda debe cumplir $c_{P}^{2} \geq 2 c_{S}^{2}$ , cosa que el constructor impone, y ese rechazo es lo primero que hay que mirar cuando un mapa de sedimento o de material compuesto no se deja construir.

3. Superficies libres y ondas de Rayleigh

Todos los lados son rígidos mientras no se declare otra cosa, y aquí rígido significa empotrado. El contorno mantiene la velocidad de partícula a cero, así que refleja con la polaridad opuesta a la de una superficie libre y no soporta onda de superficie alguna. Ese es el modelo correcto para un cuerpo que continúa en algo mucho más rígido, y el equivocado para un cuerpo en aire o en vacío, que es lo que son la mayoría de los problemas estructurales. Una placa, un bloque o una viga que se pretendan libres tienen que llevar, por tanto, todos sus lados expuestos declarados libres, y el síntoma de olvidarlo es silencioso: los frentes de onda siguen pareciendo verosímiles, pero el tren de Rayleigh que la figura de abajo hace evidente sencillamente no está, las reflexiones vuelven invertidas y no encaja ninguna forma cerrada de placa libre. La tercera opción es "absorbing", que termina un semiespacio para que el borde del dominio no se vea en absoluto; las tres se eligen por lado, que es exactamente cómo funciona el ejemplo del semiespacio: libre arriba y absorbente en los otros tres. Como no hay contrapartida elástica de FDTD2D.plot_geometry, vuelve a leer la especificación de contornos del objeto ElasticBoundaries antes de comprometerte con una ejecución larga.

Un lado declarado "free" se convierte en superficie libre de tracción mediante la imagen de esfuerzos (Moczo et al., Ec. 9.9): el esfuerzo normal se fija a cero en el plano de la superficie y el esfuerzo de cizalla por encima es su imagen antisimétrica. Esa única condición de contorno es la que hace a las superficies portadoras de ondas: golpea verticalmente la superficie libre de un bloque de aluminio (el problema de Lamb) y, además de los frentes de cuerpo P y S, una onda de Rayleigh rueda por la superficie a $c_{R} \approx 0, 93 c_{S}$ , la raíz de la ecuación característica exacta de Rayleigh (Cremer, Heckl y Petersson, Ec. 3.149). Los tests fijan la velocidad de Rayleigh medida a esa raíz dentro del 2 % y la onda de flexión de una tira delgada libre-libre a la dispersión de placa de Kirchhoff $c_{B} \propto ω$ en su dominio de placa delgada. Las superficies libres siguen la regla de muestreo de Rayleigh: reserva 15-20 celdas por longitud de onda ahí, pues la superficie de imagen de segundo orden es la parte más dispersiva del esquema.

Instantánea de la velocidad de partícula vertical en un bloque de aluminio de 0,6 por 0,3 metros unas centésimas de milisegundo después de un golpe vertical sobre su superficie libre superior: el frente P de compresión es el arco exterior, el frente S de cizalla el arco interior a aproximadamente la mitad de su radio, y los lóbulos más intensos se pegan a la superficie justo detrás del frente S, etiquetados como la onda de Rayleigh; arcos punteados marcan los radios exactos de P y S

Una instantánea, tres velocidades: el frente P ha recorrido el doble de distancia que el frente S (arcos punteados en los radios exactos $c_{P} t$ y $c_{S} t$ ), y el movimiento más intenso viaja a lo largo de la superficie libre como la onda de Rayleigh, un pelo más lenta que S.

La esponja es la herramienta equivocada allí donde vive la física de esta página. Aquí no hay una capa elástica perfectamente adaptada, y la esponja absorbe una onda de Rayleigh o de Scholte rasante mucho peor que una onda de cuerpo, precisamente porque una onda guiada por la superficie corre a lo largo de la capa en vez de meterse en ella. Tres consecuencias sobre las que conviene actuar. Un lado no puede ser libre y absorbente a la vez, así que donde una superficie libre se encuentra con una esponja en una esquina los esfuerzos imagen y la rampa de absorción trabajan uno contra otro: mantén la fuente al menos a una longitud de onda de esa esquina y termina la ejecución antes de que la onda de superficie llegue a ella. Una onda de interfase o de superficie necesita una capa más gruesa que una onda de cuerpo: empieza por aproximadamente el doble del ancho que usarías para P y S, y verifícalo alargando la ejecución y vigilando si vuelve algo. Y cuando solo importan las primeras llegadas, unos lados rígidos sin más con una ejecución más corta que el primer eco de pared salen más limpios que cualquier esponja: es exactamente lo que hace la instantánea de arriba, a 73 µs sobre un bloque de 0,6 m cuyo retorno P de la pared lateral más cercana llega a 95 µs y cuya onda de Rayleigh alcanza la esquina a 103 µs.

Mostrar el código de esta figura

import numpy as np
import matplotlib.pyplot as plt
from phonometry import simulation

# Un bloque de aluminio de 0,6 x 0,3 m, golpeado en su superficie libre.
dx, w = 0.001, 8e-6
dt = 0.6 * dx / (6320.0 * np.sqrt(2.0))
steps = round(7.3e-5 / dt)
res = simulation.elastic_fdtd_simulation(
    6320.0, 3130.0, dx, 7.3e-5, rho=2700.0, shape=(300, 600),
    sources=[simulation.ForceSource(
        ix=300, iy=0, direction="y", amplitude=1e6,
        waveform=simulation.GaussianPulse(0, 0, width=w).value)],
    boundaries=simulation.ElasticBoundaries({"top": "free"}),
    recording=simulation.ElasticRecording(snapshot_every=steps,
                                          snapshot_field="vy"),
)
res.plot(kind="snapshot", language="es")
plt.show()

La instantánea se detiene a 73 µs porque unos lados rígidos y una ejecución corta eran la manera más barata de conseguir un fotograma limpio, y un fotograma basta para las dos velocidades de cuerpo: los arcos punteados están dibujados en los radios exactos $c_{P} t$ y $c_{S} t$ , así que son una medición y no una etiqueta. Lo que un solo fotograma no puede zanjar es la tercera onda. Para ella no se dibuja nada analítico, su velocidad se separa un 7 % de la del frente de cizalla, y su propiedad definitoria, que se queda en la superficie mientras los frentes de cuerpo se marchan, es una afirmación sobre lo que pasa después.

Así que ejecuta el mismo golpe durante 163 µs, con esponjas en los otros tres lados y el dominio ensanchado hasta que la esquina donde una esponja se encuentra con la superficie libre quede fuera del encuadre, y ejecútalo dos veces: una con el lado superior declarado "free" y otra dejándolo en el valor por defecto del esquema, la pared rígida empotrada. El mismo aluminio, el mismo golpe, los mismos dos frentes de cuerpo, y ningún tren de superficie en la segunda. Eso es lo que afirma esta sección, dibujado: la onda de Rayleigh la hace la condición de contorno, no el material.

Descargar la animación (WebM)

Fíjate en lo que hacen los frentes al final. Los dos arcos de cuerpo se salen del encuadre y siguen su camino; los lóbulos de superficie no van a ninguna parte salvo a lo largo de la superficie, a su propia velocidad más lenta, y siguen ahí cuando las ondas de cuerpo ya se han marchado. El recuadro traza dos sondas de superficie separadas 0,15 m, y el retardo entre sus picos es la medición: el $c_{R}$ simulado frente a 2921 m/s, la raíz de la ecuación característica exacta de Rayleigh citada arriba, que los tests fijan dentro del 2 %.

Los dos paneles de campo están dibujados con escalas de color propias, y el clip así lo indica, porque una fuerza vertical aplicada sobre una superficie empotrada casi se cancela contra su propia imagen y esa ejecución radia mucho menos que la libre. Las trazas de sonda de debajo conservan una única escala compartida, que es donde la comparación sigue siendo cuantitativa: la llegada de superficie que falta es un pico que falta, no un pico reescalado.

La otra onda que una superficie libre hace posible, la onda de flexión de una tira delgada libre-libre, la que arriba queda fijada a la dispersión de Kirchhoff, también merece verse en movimiento. El clip de abajo es este mismo esquema lanzando un paquete de flexión de 4 kHz por una placa de acero de 10 mm: en el panel de control la placa sigue recta y el paquete simplemente se marcha, y en el panel de la unión una placa perpendicular del mismo espesor convierte la esquina en un dispersor. Allí el paquete se parte en las ondas de flexión reflejada y transmitida que la forma cerrada tasa en $τ_{12} (0°) = 0, 5$ , más el rápido precursor de onda en el plano que corre por delante por la placa receptora: la conversión de modo que el modelo de unión articulada deja fuera a propósito, y la razón por la que esta página necesita un modelo elástico y no uno de flexión. La guía de transmisión de onda de flexión desmonta esa misma ejecución frente al índice de reducción vibratoria de EN 12354.

Descargar la animación (WebM)

4. Acoplamiento fluido-sólido a incidencia normal

Como el fluido es simplemente el límite c_s = 0 del esquema elástico, un contacto agua-acero es una sola simulación con dos bandas en los mapas de material, y toda la física del acoplamiento fluido-sólido emerge de la misma actualización. Las constantes Material y ElasticFDTD2D.from_regions construyen esos mapas por capas sin manipular los arrays a mano: un material de fondo y una lista de regiones pintadas.

Un Material son tres números medibles (la velocidad de compresión, la velocidad de cizalla y la densidad), y una velocidad de cizalla nula marca un fluido. Se distribuyen cinco constantes, y de ellas sale cualquier otro número de esta página:

Constante	$c_{P}$ [m/s]	$c_{S}$ [m/s]	$ρ$ [kg/m³]
`AIR`	343	0	1,2
`WATER`	1480	0	1000
`CONCRETE`	3800	2250	2400
`STEEL`	5900	3200	7850
`ALUMINIUM`	6320	3130	2700

Son cifras redondas nominales, para documentación y validación, no datos certificados, así que una probeta real hay que introducirla a mano: por eso la sección 1 escribe los 6320 / 3130 / 2700 del bloque de aluminio en vez de importar ALUMINIUM. El constructor impone una sola cota, $c_{P}^{2} \geq 2 c_{S}^{2}$ , que es la afirmación de que el primer parámetro de Lamé no es negativo o, equivalentemente, de que el coeficiente de Poisson no es negativo. Un trío rechazado suele significar una velocidad de cizalla introducida demasiado alta, o una velocidad de compresión medida sobre una barra delgada en vez de en el seno del material: la de barra delgada es la menor de las dos, y sustituir una por otra es un error frecuente.

import numpy as np
from phonometry.simulation import ElasticFDTD2D, STEEL, WATER

# Una columna de agua de 12 m sobre un semiespacio de acero, como tira 1D.
dx = 0.005
sim = ElasticFDTD2D.from_regions(
    (2400, 3), dx, background=WATER, regions=[(np.s_[1200:, :], STEEL)])
y = (np.arange(2400) + 0.5) * dx
p0 = np.exp(-(((y - 3.0) / 0.15) ** 2))[:, None]     # pulso plano
sim.txx[:] = -p0
sim.tyy[:] = -p0
y_face = np.arange(1, 2400) * dx + 0.5 * WATER.c_p * sim.dt
sim.vy += (np.exp(-(((y_face - 3.0) / 0.15) ** 2))
           / (WATER.rho * WATER.c_p))[:, None]       # unidireccional, consistente con el leapfrog
trace = []
for _ in range(round(3.4e-3 / sim.dt)):
    sim.step()
    trace.append(sim.p[900, 1])
trace = np.asarray(trace)
t = (np.arange(trace.size) + 1) * sim.dt
incident = trace[t < 1.8e-3].max()
echo = trace[t > 2.6e-3]
print(round(float(echo[np.abs(echo).argmax()] / incident), 3))   # 0.938
# (Z2 - Z1)/(Z2 + Z1) = 0.938: a incidencia normal no se excita cizalla
# y el acero se comporta como un líquido de su rho y c_P (B&G Ec. 4.2.27)

Ese 0,938 es una razón entre dos eventos de una misma traza, y la traza es donde se ve si la ejecución fue lo bastante limpia como para citarlo. ElasticFDTDResult.plot() dibuja exactamente esto por defecto: historias de sonda, no mapas de campo; las instantáneas de las secciones 3 y 6 son .plot(kind="snapshot").

Historia de presión en una sonda dentro de una columna de agua de 12 metros sobre un semiespacio de acero, a lo largo de 3,4 milisegundos: un único pulso gaussiano limpio de amplitud unidad que pasa por la sonda a 1,02 milisegundos camino abajo, un tramo plano y silencioso, y un segundo pulso de la misma polaridad y 0,938 veces la amplitud a 3,04 milisegundos, la reflexión que vuelve del acero, con los dos instantes de llegada marcados y la razón anotada frente a la fórmula de impedancias en forma cerrada

Toda la medición en un solo eje. La sonda está 7,5 m por debajo de la fuente y 7,5 m por encima del contacto, así que el pulso incidente pasa por ella a 1,02 ms y la reflexión vuelve a 3,04 ms: dos instantes fijados por la geometría y por los 1480 m/s del agua, y los dos merecen comprobarse antes de leer las amplitudes. Los dos eventos conservan la misma polaridad, que es lo que hace una reflexión en un medio más duro, y su razón es 0,938 frente al valor exacto $(Z_{2} - Z_{1}) / (Z_{2} + Z_{1}) = 0, 938$ . El tramo plano entre ambos es la prueba de que no está contribuyendo nada más del dominio.

Mostrar el código de esta figura

import matplotlib.pyplot as plt

# La ejecución de arriba, releída: `res.plot()` sobre un resultado congelado
# dibuja esto directamente. Aquí la traza se acumuló a mano, así que se
# dibuja a mano.
plt.plot(t * 1e3, trace)
plt.axvline((4.5025 - 3.0) / WATER.c_p * 1e3, ls=":")   # incidente, 1,02 ms
plt.axvline((3.0 + 1.4975) / WATER.c_p * 1e3, ls=":")   # eco, 3,04 ms
plt.xlabel("Tiempo [ms]")
plt.ylabel("Presión [Pa]")
plt.show()

5. Incidencia oblicua: conversión de modo

Fuera de la incidencia normal el sólido deja de ser un líquido. El sonido incidente se refracta en dos ondas transmitidas, P y SV (conversión de modo), y la refracción obedece la ley de Snell por separado para cada una de ellas, así que cada una tiene su propio ángulo crítico: $arcsin (c_{1} / c_{P 2})$ para la onda de compresión y $arcsin (c_{1} / c_{S 2})$ para la de cizalla. Ambos salen solo de las velocidades de onda (la densidad entra únicamente en las amplitudes), así que son dos líneas de aritmética sobre las constantes de material de la sección 4.

import numpy as np

# Los dos ángulos críticos salen solo de las velocidades. Sustituye tu
# propio par de Material y esas mismas dos líneas te lo responden.
print(round(float(np.degrees(np.arcsin(WATER.c_p / STEEL.c_p))), 1))   # 14.5
print(round(float(np.degrees(np.arcsin(WATER.c_p / STEEL.c_s))), 1))   # 27.5

El dibujo de abajo reúne las tres disposiciones que describen esta sección y las dos siguientes, sobre una misma línea de interfase y a tres incidencias, junto con aquello que un error de una fila estropea: dónde está físicamente el contacto sobre la malla escalonada.

Tres paneles que comparten una misma línea horizontal de interfase fluido-sólido, con agua arriba y sólido rayado abajo. Panel a, incidencia normal: una flecha incidente hacia abajo, una flecha reflejada hacia arriba y una única flecha P transmitida, con el coeficiente de reflexión 0,938 indicado debajo. Panel b, incidencia oblicua: un rayo incidente formando ángulo con la normal, su rayo reflejado y dos rayos transmitidos etiquetados P y SV que salen con ángulos distintos, con una sonda dibujada sobre la normal 0,12 metros por encima del punto de impacto y los ángulos críticos de 14,5 y 27,5 grados indicados debajo. Panel c, el caso de Scholte: una explosión diez metros por encima de un contacto agua-sedimento y un tren oscilante pegado a la interfase a 1436 metros por segundo, con envolventes evanescentes discontinuas a ambos lados. Un detalle ampliado debajo muestra la última fila de fluido, la primera fila de sólido y el plano de caras entre ambas.

Esos dos ángulos parten el rango en cuatro regímenes. La columna titulada fluido equivalente es esa misma interfase modelada como un fondo sin cizalla del mismo $c_{P}$ y $ρ$ , el modelo que lleva el módulo submarino, y la diferencia entre las dos columnas es toda la razón de ser de este esquema.

Incidencia	Ondas transmitidas	\|V\| elástico	\|V\| fluido equivalente	Qué significa
0°	solo P	0,938	0,938	no se excita cizalla, así que el acero es exactamente un líquido de su $ρ$ y su $c_{P}$
de 0° a 14,5°	P y SV	0,938, que baja a 0,937 en 9° y luego sube hasta el ángulo crítico	0,938, que sube hasta 1	la cizalla ya se lleva su parte de la potencia transmitida, y el modelo fluido ya devuelve de más
de 14,5° a 27,5°	solo SV, P evanescente	baja hasta 0,905 en 16,4°; 0,918 en 20°	1 (reflexión total)	el modelo fluido sobrestima la energía devuelta en toda la fuga por cizalla
más allá de 27,5°	ninguna se propaga	1, con una fase	1, con una fase distinta	reflexión total en ambos casos; la diferencia de fase es aquello a lo que se engancha una onda de Scholte

Lee dos veces la tercera fila: entre los ángulos críticos un fluido equivalente lo refleja todo, mientras que la interfase real deja escapar cerca de una décima parte de la amplitud en forma de onda de cizalla (|V| baja hasta 0,905 en 16,4°). La batería de validación lanza haces oblicuos con portadora a 10° y 20° y reproduce el coeficiente de reflexión exacto de Brekhovskikh y Godin (Ecs. 4.2.22-4.2.26) dentro de una fracción de punto porcentual, conversión de modo incluida, así que la columna elástica de arriba está medida y no solo calculada.

Cómo se hace esa medición, porque es una técnica que merece la pena robar. Un haz inclinado con portadora gaussiana se escribe como condición inicial sobre txx y tyy, con la velocidad de partícula medio paso temporal atrás, consistente con el leapfrog, repartida entre vx y vy por $sin θ$ y $cos θ$ : el mismo truco unidireccional que la sección 4 usa a incidencia normal, girado. La escena se ejecuta después dos veces, una con el sólido pintado y otra con el fluido continuado de largo a través de él, y la ejecución de referencia se fuerza al mismo paso temporal escalando su número de Courant por $c_{agua} / c_{acero}$ , de modo que las dos trazas se restan muestra a muestra. La diferencia es el campo reflejado a solas, artefactos de pared lateral incluidos: el análogo elástico de la técnica de sustracción de ISO 13472-1. La sonda está justo encima del punto donde el eje del haz corta la interfase, 0,12 m por encima, de manera que los centros de los haces incidente y reflejado la cruzan con desplazamientos transversales simétricos y la envolvente se cancela en la razón; $∣ V ∣$ es entonces la razón espectral ponderada por amplitud sobre la banda de la portadora, no una lectura de pico a pico.

El dibujo de arriba es un diagrama de rayos: dice adónde va cada onda. Lo que un rayo no puede decir es qué aspecto tiene una «P evanescente», y en eso consiste entera la tercera fila de la tabla. El clip de abajo ejecuta la escena tres veces, con una incidencia de cada uno de los tres primeros regímenes (10°, por debajo de ambos ángulos críticos; 20°, entre ellos; 35°, pasado el de cizalla), sobre el mismo semiespacio de WATER y STEEL, lanzando contra el contacto un haz sostenido con gradiente de fase. Léelo como un retrato de los regímenes, no como una segunda medición: cada panel lleva anotado el $∣ V ∣$ en forma cerrada de la tabla de arriba, porque el haz que cabe en un panel de 90 mm mide seis longitudes de onda de ancho, y con ese ancho la sonda queda en el campo próximo de la fuente en lugar de en un haz ya formado. Ensancharlo hasta que la lectura fuese de fiar exigiría un panel cuatro veces más ancho, en el que las franjas de interferencia serían demasiado finas para verse. La columna medida de la tabla sale, en cambio, de la batería de validación, que sí tiene sitio para hacerlo como es debido.

Descargar la animación (WebM)

A 10° salen del contacto dos haces por las direcciones de Snell dibujadas a trazos, y salen con ángulos distintos porque viajan a velocidades distintas: la onda de cizalla, más lenta, se hunde a 22,1° de la normal, y la de compresión, más rápida, se abre mucho más, hasta 43,8°. A 20° el haz P no es un haz más rasante: ha desaparecido. Lo que ocupa su lugar es una piel pegada al contacto que muere en menos de una longitud de onda y no lleva nada hacia abajo, mientras el haz de cizalla se lleva él solo lo que cruza; y esa fuga por cizalla es la razón de que la reflexión se quede en 0,918 en lugar del 1,000 que devolvería un fluido equivalente. A 35° los dos son evanescentes y nada se propaga alejándose por el acero: lo que llena el sólido ahí es una onda inhomogénea que corre a lo largo del contacto y muere con la profundidad, brillante y sin ir a ninguna parte. Encima del contacto, mientras tanto, el haz incidente y una reflexión total en amplitud pero desplazada en fase se resuelven en un patrón estacionario. De esa fase construye la sección siguiente una onda de Scholte.

Cada panel dibuja su mitad de acero con su propia ganancia de representación, anotada, y la razón merece conocerse: el agua excita el acero a través de un salto de impedancia de treinta veces, así que el campo transmitido lleva una velocidad de partícula pequeña incluso allí donde lleva potencia real, y los tres paneles se diferencian entre sí de un modo que ninguna ganancia única abarca. A 10° la energía transmitida se marcha en dos haces, así que la amplitud es pequeña en cualquier punto; a 35° no se marcha nada y el campo evanescente se amontona contra el contacto, brillante y sin ir a ninguna parte. El brillo en el acero, por tanto, no es potencia. Para eso, lee el $∣ V ∣$ impreso y el factor de ganancia de cada panel.

Esta es exactamente la física de cizalla que el reflection_coefficient de fondo fluido del módulo submarino no puede representar, y la razón por la que un lecho marino con cizalla pierde más energía que su modelo de fluido equivalente.

6. La onda de interfase de Scholte

La reflexión total más allá del ángulo crítico S deja ambos medios con campos solo evanescentes, y esos campos pueden acoplarse en una verdadera onda de interfase: la onda de Scholte, el análogo fluido-sólido de la onda de Rayleigh. Viaja más despacio que la velocidad del sonido en el fluido y que la de cizalla del sólido, no tiene corte en baja frecuencia y no es dispersiva sobre semiespacios homogéneos; scholte_speed resuelve su ecuación característica exacta (B&G Ec. 4.4.20).

from phonometry.simulation import Material, STEEL, WATER, scholte_speed

print(round(scholte_speed(WATER, STEEL), 1))       # 1479.6, un 0,03 % por debajo de la del agua
seabed = Material(c_p=3500.0, c_s=2000.0, rho=2500.0)
print(round(scholte_speed(Material(1500.0, 0.0, 1000.0), seabed), 1))  # 1436.0

Los dos casos usan aguas distintas a propósito. WATER es la constante de la biblioteca, a 1480 m/s, mientras que el caso del sedimento sigue el benchmark de van Vossen, cuyo fluido va a 1500 m/s: así que los dos déficits de abajo, el 0,03 % y el 4 %, están medidos cada uno contra su propia agua, y 1436 m/s no debe compararse nunca con 1480 m/s.

Esos dos números cuentan toda la historia. Sobre un lecho rígido (acero) la onda de Scholte se pega a la velocidad del agua (1480 m/s) hasta un 0,03 % y su cola evanescente alcanza ~7 longitudes de onda dentro del agua: es en esencia una onda de agua rasante y no puede separarse por tiempo de vuelo en ningún dominio razonable (sobre contactos aire-sólido el déficit se desploma a ~ $1 0^{- 12}$ de $c$ y la onda es directamente inobservable). Sobre un sedimento blando la velocidad cae muy por debajo de la del agua y la onda queda confinada a menos de media longitud de onda del contacto, y por eso las ondas de interfase medidas en el lecho marino son una herramienta estándar para sondear la velocidad de cizalla del sedimento. La instantánea de abajo ejecuta ese caso blando: una explosión justo sobre el fondo, y el rasgo tardío más intenso es el tren de Scholte pegado a la interfase, cronometrado por la batería de tests a 1436 m/s dentro del 2 % entre dos sondas en el contacto.

Instantánea de la velocidad de partícula vertical en un dominio de 500 por 200 metros de agua sobre un semiespacio de sedimento blando, 0,2 segundos después de una explosión 10 metros por encima del contacto: el frente directo del agua se arquea hacia arriba y a la derecha, frentes oblicuos radian hacia el sedimento, y los lóbulos más intensos forman un tren compacto pegado a la línea punteada de la interfase, etiquetado como la onda de Scholte, evanescente a ambos lados

Los medios del benchmark de van Vossen (agua sobre un sedimento 3500/2000/2500): el tren de Scholte repta por el contacto a 1436 m/s, evanescente hacia ambos medios, mientras la onda directa del agua corre por delante a 1500 m/s.

Mostrar el código de esta figura

import numpy as np
import matplotlib.pyplot as plt
from phonometry import simulation

# Agua sobre un lecho blando; Ricker explosivo de 50 Hz a 10 m del contacto.
ny, nx, dx = 200, 500, 1.0
c_p = np.full((ny, nx), 1500.0)
c_s = np.zeros((ny, nx))
rho = np.full((ny, nx), 1000.0)
c_p[100:], c_s[100:], rho[100:] = 3500.0, 2000.0, 2500.0
f0, t0 = 50.0, 0.030
dt = 0.6 * dx / (3500.0 * np.sqrt(2.0))
steps = round(0.232 / dt)

def ricker(t):
    a = (np.pi * f0 * (t - t0)) ** 2
    return (1.0 - 2.0 * a) * np.exp(-a)

res = simulation.elastic_fdtd_simulation(
    c_p, c_s, dx, 0.232, rho=rho,
    sources=[simulation.ExplosionSource(ix=60, iy=89, waveform=ricker,
                                        amplitude=1e3)],
    boundaries=simulation.ElasticBoundaries("absorbing",
                                            absorbing_layer_cells=20),
    recording=simulation.ElasticRecording(snapshot_every=steps,
                                          snapshot_field="vy"),
)
res.plot(kind="snapshot", language="es")
plt.show()

7. Placas sumergidas y la ley de masas

Una placa sumergida cierra el lazo con el módulo de acústica de la edificación. A incidencia normal no se excita cizalla, así que una placa elástica en agua es exactamente la capa fluida de tres medios de B&G §2.4: su transmisión sigue la forma cerrada de la Ec. 2.4.14, que se reduce a la conocida ley de masas para placas delgadas y bajas frecuencias y predice transmisión total en las resonancias de espesor de media onda $f_{n} = n c_{P} / (2 h)$ (Ec. 2.4.19). Para una placa de acero de 10 mm esa primera resonancia cae en 295 kHz, y una sola ejecución FDTD de banda ancha reproduce toda la curva: 5,8 dB a 10 kHz (donde mass_law_transmission_loss con agua como fluido ambiente coincide con la forma exacta hasta 0,02 dB), 18,1 dB a 50 kHz, y un mínimo de pérdida de transmisión dentro del 0,1 % de la resonancia de 295 kHz. La misma batería lleva al límite el contraste extremo de un contacto aire-acero (razón de impedancias ~ $1 0^{5}$ :1): estable durante 10 000 pasos con la amplitud reflejada conservada al 0,5 %.

Dibujo del montaje de la ejecución de transmisión con placa sumergida. A la izquierda, una tira vertical alta y estrecha dibujada a escala: 0,75 metros de agua con una placa de acero de 10 milímetros a 0,35 metros, un pulso plano unidireccional lanzado hacia abajo a 0,25 metros, la sonda A 30 milímetros por encima de la placa y la sonda B 50 milímetros por debajo, con una nota de que la tira tiene tres celdas de ancho, un problema unidimensional ejecutado en un esquema bidimensional a medio milímetro por celda. A la derecha, las dos trazas de sonda sobre un eje temporal común de 0 a 160 microsegundos: la sonda A muestra el pulso incidente a 47 microsegundos dentro de una puerta sombreada de 65 microsegundos y su reflexión en la placa a 88 microsegundos, ya fuera de ella; la sonda B muestra solo la cola transmitida. Debajo, la fórmula de la pérdida de transmisión y un esquema enmarcado del patrón estacionario de media onda dentro de la placa con su primera resonancia de 295 kilohercios.

La medición que hay detrás de esos tres números. La escena es una tira de tres celdas de ancho y 0,75 m de alto con $Δ x = 0, 5$ mm (un problema unidimensional ejecutado en un esquema bidimensional) llena de WATER, con una región de STEEL de 10 mm pintada a 0,35 m mediante from_regions. Un pulso plano gaussiano de 1,5 mm de anchura se escribe unidireccional igual que en la sección 4, centrado en 0,25 m; esa anchura es la que lleva energía aprovechable hasta unos 340 kHz, que es lo que mete la resonancia de 295 kHz dentro de la banda medida. Una sonda está 30 mm por encima de la placa y otra 50 mm por debajo de su cara inferior, y la ejecución dura 0,45 ms.

La medición en sí es una ventana temporal, y la geometría te la regala. La sonda superior ve el pulso incidente hacia los 47 µs (0,07 m de agua a 1480 m/s) y su reflexión en la placa hacia los 88 µs (0,13 m); las dos llegadas duran unos 10 µs, así que una puerta de 65 µs aísla la incidente exactamente, no aproximadamente. La sonda inferior se mantiene libre de ecos durante todo el registro, así que no necesita puerta alguna. Después, $TL (f) = 20 lo g_{10} ∣ I (f) / T (f) ∣$ a partir de los dos espectros de magnitud. Dos comprobaciones que el lector puede hacer sobre su propia placa: el valor a 10 kHz frente a mass_law_transmission_loss con agua como fluido ambiente, y el mínimo frente a $f_{n} = n c_{P} / (2 h)$ . Esa misma lógica de ventana temporal es la que aplica una medición in situ de ISO 13472-1 con una ventana de Adrienne, así que el truco se traslada tal cual a un banco real.

Para qué sirve en realidad esa cifra de 295 kHz. Está tres décadas por encima de la banda de la acústica de la edificación, así que un panel a frecuencias de audio está siempre en el límite de placa delgada, donde la transmisión de la capa se desploma sobre la ley de masas y la resonancia de espesor puede ignorarse por completo: por eso justamente las formas cerradas de la guía de aislamiento de paneles no la mencionan nunca. La resonancia se gana el sueldo en su propia banda: la transmisión total en el espesor de media onda es el principio en que se basan la medida ultrasónica de espesores y la elección de las capas de acoplamiento de los transductores, y por eso la validación lleva la placa a propósito hasta los centenares de kilohercios. La frecuencia donde se encuentran los dos regímenes es una divisoria válida para cualquier placa: $c_{P} / (2 h)$ , que vale 295 kHz para 10 mm de acero y 59 kHz para 50 mm del mismo.

A incidencia oblicua la física de la placa se enriquece, y el clip de abajo es este mismo esquema excitando la misma placa de acero de 10 mm, ahora tendida en aire, con una onda plana sostenida que llega a 45°. Los dos paneles se diferencian en un solo número: la frecuencia de excitación, $f_{c} /2 = 603$ Hz a la izquierda y $2 f_{c} = 2413$ Hz a la derecha, a un lado y otro de la frecuencia de coincidencia de 1206 Hz que la biblioteca calcula con los mismos $m^{''}$ y $B^{'}$ usados más arriba. Todo lo demás (la placa, el ángulo, la malla y la escala de color) se mantiene fijo. Por debajo de $f_{c}$ la placa refleja casi todo y el nivel transmitido cae sobre la ley de masas oblicua; por encima, la longitud de onda de traza acústica coincide con la longitud de onda de flexión libre, la placa rerradia un haz a 45° que crece a lo largo del tramo iluminado y el nivel transmitido se mantiene en la cifra de baja frecuencia, justo donde la ley de masas exigía 12 dB más de bloqueo. El aire bajo la placa se dibuja en los dos paneles con la ganancia de representación que _weak_field_gain mide sobre el campo asentado de las dos ejecuciones juntas (×150, es decir, +44 dB) y que se imprime en el lienzo: lee las anotaciones para los niveles, no el brillo. La guía de predicción del aislamiento de paneles desmonta esa misma ejecución frente al método de la meseta y a la ley de masas.

Descargar la animación (WebM)

Qué cubre esta guía

Cubierto
El compañero elástico P-SV del método FDTD acústico: el esquema velocidad-esfuerzo de Virieux (1986) con los parámetros efectivos y la superficie libre por imagen de esfuerzos de Moczo et al. (2007), expuesto a través de elastic_fdtd_simulation/ElasticFDTD2D con ExplosionSource y ForceSource, validado contra las velocidades de las ondas de cuerpo, la raíz exacta de Rayleigh, la dispersión de flexión de Kirchhoff y el límite acústico bit a bit. El acoplamiento fluido-sólido de las secciones 4-7, con las constantes Material, ElasticFDTD2D.from_regions y scholte_speed: la reflexión oblicua con conversión de modo contra Brekhovskikh y Godin (Ecs. 4.2.22-4.2.26) en ángulos subcríticos y entre críticos, la velocidad de la onda de Scholte contra la ecuación característica exacta (Ec. 4.4.20) y por tiempo de vuelo en el benchmark de fondo blando de van Vossen (2002), la transmisión de la placa sumergida con su primera resonancia de espesor contra la forma cerrada de tres medios (Ecs. 2.4.14-2.4.19), la ley de masas del panel sumergido y la prueba de robustez con el contraste extremo aire-acero.
No cubierto
Todo lo que el esquema acústico 2D no puede hacer, este tampoco: el dominio es una sección transversal (propagación cilíndrica $1/ r$ , no $1/ r$ ), así que los niveles absolutos y las tasas de decaimiento no son los de un problema 3D. El sólido es isótropo y puramente elástico: sin anisotropía ni amortiguamiento viscoelástico más allá de la tasa de decaimiento volumétrica, así que los factores de pérdidas del material no se modelan. No hay una capa elástica perfectamente adaptada: la esponja absorbe las ondas de Rayleigh rasantes peor que las ondas de cuerpo. La dispersión de ondas guiadas más allá del régimen de placa delgada de Kirchhoff (modos de Lamb completos) es observable en los campos pero no tiene aquí oráculo dedicado en forma cerrada.

Véase también

Simulación de ondas FDTD 2D: el esquema acústico que esta guía extiende, con la malla escalonada, la cota de Courant y la maquinaria de fuentes, sondas y contornos que ambos comparten.
Predicción del aislamiento de paneles: las formas cerradas de la ley de masas y la coincidencia; la ejecución de la placa sumergida del §7 es su contraste de onda completa.
Transmisión de onda de flexión en uniones de placas: las formas cerradas de ondas de flexión del modelo de edificación; la dispersión de Kirchhoff validada en el §3 es la onda que transportan.
Propagación submarina del sonido: el modelo de reflexión de Rayleigh de fondo fluido cuya física de cizalla ausente cuantifica el §5.
Referencia de la API: simulation.elastic-fdtd.

Respuestas rápidas

¿Cuándo necesito el esquema elástico en lugar del acústico?

Siempre que alguna parte del dominio soporte cizalla: la conversión de modo a incidencia oblicua, las ondas de Rayleigh en superficies libres, las ondas de Scholte en contactos fluido-sólido y la transmisión de placas necesitan las ecuaciones elásticas. Si todas las celdas son fluidas, usa el esquema acústico: con un mapa uniforme de c_s = 0 el elástico lo reproduce bit a bit, a aproximadamente el doble de coste en memoria y en avance temporal.

¿Qué resolución de malla necesita una simulación elástica?

Resuelve al menos 10 celdas por longitud de onda más corta con la velocidad de onda más lenta del dominio, la de cizalla no nula más lenta o una velocidad del sonido de fluido aún más lenta si la hay: $Δ x \leq c_{min} / (10 f)$ . Reserva 15-20 celdas por longitud de onda a lo largo de superficies libres (el contorno de imagen de esfuerzos es la parte más dispersiva del esquema), y al menos 15 puntos por longitud de onda cuando haya que cronometrar una onda de interfase (van Vossen et al. 2002).

¿Qué es una onda de Scholte y a qué velocidad viaja?

La onda de Scholte es la verdadera onda de interfase de un contacto fluido-sólido: evanescente hacia ambos medios, con movimiento elíptico de partícula, sin corte en baja frecuencia y no dispersiva sobre semiespacios homogéneos. Su velocidad resuelve la ecuación característica exacta (Brekhovskikh y Godin, Ec. 4.4.20, scholte_speed) y siempre queda por debajo de la velocidad del sonido en el fluido y de la de cizalla del sólido: un 0,03 % por debajo de la del agua (1480 m/s) sobre acero (1479,6 m/s), pero un 4 % por debajo de los 1500 m/s del agua oceánica del ejemplo de sedimento (1436 m/s para esa agua sobre un lecho 3500/2000/2500), y por eso las ondas de interfase del lecho marino sondean la rigidez a cizalla del sedimento.

Referencias

Brekhovskikh, L. M. y Godin, O. A. (1990). Acoustics of layered media I: Plane and quasi-plane waves. Springer. https://doi.org/10.1007/978-3-642-52369-4Los oráculos fluido-sólido de las secciones 4-7: el coeficiente de reflexión oblicua con conversión de modo (Ecs. 4.2.22-4.2.26) y sus límites en los ángulos críticos (Ecs. 4.2.27-4.2.31), la ecuación característica exacta de Scholte (Ecs. 4.4.18-4.4.20) con las asintóticas de fluido ligero (Ecs. 4.4.21-4.4.24), y la transmisión de capa de tres medios (Ecs. 2.4.10-2.4.19).
Cremer, L., Heckl, M. y Petersson, B. A. T. (2005). Structure-borne sound (3.ª ed.). Springer. https://doi.org/10.1007/b137728Los oráculos analíticos del esquema elástico: la ecuación característica exacta de Rayleigh (Ec. 3.149) y la dispersión de flexión de Kirchhoff con su corrección de espesor (Ecs. 3.83-3.89, 3.196b).
Moczo, P., Kristek, J., Galis, M., Pazak, P. y Balazovjech, M. (2007). The finite-difference and finite-element modeling of seismic wave propagation and earthquake motion. Acta Physica Slovaca, 57(2), 177-406. Los parámetros efectivos de malla en medio heterogéneo de la sección 2 (módulo de cizalla armónico, densidad aritmética; Ecs. 7.37-7.39) y la superficie libre por imagen de esfuerzos de la sección 3 (Ec. 9.9).
van Vossen, R., Robertsson, J. O. A. y Chapman, C. H. (2002). Finite-difference modeling of wave propagation in a fluid-solid configuration. Geophysics, 67(2), 618-624. https://doi.org/10.1190/1.1468623El benchmark fluido-sólido del mismo esquema escalonado: las medias armónica del módulo de cizalla y aritmética de la densidad (Ecs. 9-10) que satisfacen implícitamente las condiciones de interfase, la configuración de fondo blando del caso de Scholte de la sección 6 y la regla de >= 15 puntos por longitud de onda para el esquema O(2,2) con una onda de interfase.
Virieux, J. (1986). P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics, 51(4), 889-901. https://doi.org/10.1190/1.1442147El esquema elástico de la sección 1: el sistema velocidad-esfuerzo (Ec. 2), la celda totalmente escalonada y su actualización (Fig. 1, Ec. 5), la cota de Courant solo con P (Ecs. 6-7), las relaciones de dispersión (Ecs. 13-14) y el líquido como límite sin cizalla.

Creado y mantenido por José Manuel Requena Plens ORCID Todos los proyectos