Source: https://jmrplens.github.io/phonometry/es/simulation/fdtd-simulation/
La mayor parte de phonometry predice un **número**: un nivel, un tiempo de
reverberación, una pérdida por transmisión. El dominio `simulation` calcula el
**propio campo de ondas**: `fdtd_simulation` integra las ecuaciones acústicas
lineales sobre una malla 2D con el método de **diferencias finitas en el
dominio del tiempo (FDTD)**, de modo que la reflexión, la difracción, la
interferencia, la refracción en medios inhomogéneos y el comportamiento modal
emergen de primeros principios en lugar de modelarse término a término. La
implementación sigue la formulación de referencia para sonido en exteriores de
Attenborough y Van Renterghem, *Predicting Outdoor Sound* (2.ª ed., CRC Press
2021), capítulo 4: el esquema presión-velocidad escalonado en espacio y en
tiempo (Ecs. 4.11-4.12), la condición de estabilidad de Courant
(Ecs. 4.13-4.14), los contornos rígidos como velocidad normal nula en la cara
(Ec. 4.32) y el contorno de impedancia real independiente de la frecuencia
(Ecs. 4.33-4.35).
El esquema es **determinista por diseño**: aritmética float64, sin
números aleatorios y avance numpy monohilo, de modo que las mismas entradas
producen salidas idénticas bit a bit en la misma plataforma. Es el motor de
las animaciones FDTD de esta documentación, elevado a API pública con
fuentes, sondas de presión, obstáculos rasterizados, condiciones de contorno
por lado y un objeto de resultado congelado.
Aquí tienes una de esas animaciones, y es un anuncio honesto de lo que
construye el resto de la página. En ella no hay nada dibujado: la columnata es
una `obstacle_mask` booleana de círculos rasterizados y el frente de onda es un
único paquete de onda plana unidireccional con una envolvente gaussiana de una
longitud de onda de ancho, lanzado en $x = 0{,}30$ m dentro de una sala de
4 m × 1 m de paredes rígidas cuyos dos extremos absorben mediante esponjas
ocultas fuera del encuadre. La portadora es de 800 Hz, así que la longitud de
onda es de 42,9 cm y las columnas de 10 a 17 cm son aproximadamente de un cuarto
a dos quintos de ella: el régimen en el que un cilindro rígido a la vez proyecta
una sombra legible y rerradia con fuerza, que es la razón por la que la cola que
llena la sala está estructurada y no es ruido. Esa cola es dispersión múltiple
determinista: es el aspecto que tiene un campo difuso *antes* de hacer ninguna
hipótesis estadística sobre él. La malla es el ejemplo resuelto de la regla que
deduce la [sección 6](#dispersion-numerica-y-precision): el hueco más
estrecho de esta disposición es de 6,6 cm entre una columna y una pared, así que
$\Delta x = \min(\text{apertura más pequeña}/4,\ \lambda/8)$ permite hasta
1,6 cm, y el clip corre a 2,5 mm porque un banner necesita la definición.
## 1. El esquema: una ecuación de onda sobre una malla
En un medio en reposo, las ecuaciones linealizadas de la dinámica de fluidos
se reducen a un sistema de primer orden en la presión acústica $p$ y la
velocidad de partícula $\mathbf{v}$ (Attenborough y Van Renterghem,
Ecs. 4.3-4.4):
$$
\frac{\partial p}{\partial t} + \rho c^2\,\nabla\!\cdot\!\mathbf{v} = 0,
\qquad
\frac{\partial \mathbf{v}}{\partial t} + \frac{1}{\rho}\,\nabla p = 0 .
$$
FDTD discretiza ambas sobre una **malla escalonada** (el análogo acústico de
la celda de Yee): la presión vive en los centros de celda y cada componente
de velocidad en las caras, a media celda, y los dos campos avanzan en
**leapfrog**, medio paso temporal desfasados (Ecs. 4.11-4.12). Evaluar cada
gradiente espacial exactamente donde lo necesita el otro campo cuadruplica la
precisión frente a una malla colocalizada (Ec. 4.9 frente a 4.10) y permite
actualizaciones in situ. Como solo se almacenan las caras interiores, el
borde del dominio es una **pared perfectamente rígida** (velocidad normal
nula, Ec. 4.32) salvo que se pida otro contorno.
El esquema explícito solo es estable mientras un frente de onda cruce como
mucho una celda por paso temporal. Con celdas cuadradas, el **número de
Courant** (Ec. 4.13) es
$$
\mathrm{CN} = c\,\Delta t\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}
= \frac{c\,\Delta t\,\sqrt{2}}{\Delta x} \le 1,
$$
y `fdtd_simulation` deriva el paso temporal del parámetro `cfl` (este
$\mathrm{CN}$, por defecto 0,6) y de la mayor velocidad del sonido del mapa. El
esquema solo es **condicionalmente estable**: con $\mathrm{CN} \le 1$ la
actualización ni crea ni destruye energía, y por encima del límite todo modo
de longitud de onda próxima a la escala de la malla crece exponencialmente,
así que la ejecución llega a `inf` en unos pocos centenares de pasos en vez
de degradarse con suavidad (Ec. 4.14). Por eso `fdtd_simulation` rechaza
cualquier `cfl` fuera de $(0, 1)$ en lugar de dejar que la ejecución
arranque. La sección 6 usa además un número de Courant
*por eje*, $S = c\,\Delta t/\Delta x = \mathrm{CN}/\sqrt{2}$, en la relación de
dispersión; los dos no son el mismo número, y el `cfl = 0.6` por defecto
significa $S = 0{,}424$.
```python
from phonometry import simulation
# Un dominio de aire de 3.0 x 2.0 m: 300 x 200 celdas de 1 cm.
res = simulation.fdtd_simulation(
343.0, # c: velocidad del sonido [m/s], escalar o mapa sobre la malla
0.01, # dx: tamaño de celda cuadrada [m]
2.0e-3, # duration: tiempo simulado [s] (dt se deriva, no se da)
shape=(200, 300), # (ny, nx) celdas -> un dominio de 2.0 x 3.0 m
sources=[simulation.GaussianPulse(ix=60, iy=100, width=3.0e-4)],
probes=[(200, 100)],
)
print(res.size) # (3.0, 2.0) metros
print(round(res.dt * 1e6, 2)) # 12.37 microsegundos (CN = 0.6)
res.plot() # historias de presión en las sondas (figura de la sección 5)
```
El paso temporal no es una entrada. Sale de `cfl` y de la mayor velocidad
del sonido del mapa, y `res.dt` informa del valor que se usó: por eso la
duración se da en segundos de tiempo simulado y nunca en pasos.
La malla es de índices: la celda ``(ix, iy)`` tiene su centro en
``((ix + 0.5) * dx, (iy + 0.5) * dx)`` metros, con las filas dibujadas hacia
abajo (la convención de ``imshow``), de modo que una posición en metros se
convierte con ``ix = round(x / dx - 0.5)``.
## 2. Fuentes, sondas, obstáculos y contornos
Tres tipos de fuente inyectan una **fuente blanda** (una contribución de
presión aditiva que no dispersa las ondas que pasan) en una celda de la
malla: `GaussianPulse` (un pulso de banda ancha con semianchura temporal
`width`), `CWSource` (un tono sinusoidal con rampa de coseno alzado para que su
arranque no salpique el campo con un transitorio de banda ancha) y
`SignalSource` (una forma de onda muestreada arbitraria, interpolada
linealmente sobre los pasos temporales de la simulación). Las sondas
registran la presión en su celda en cada paso temporal dentro del resultado.
**Ondas planas «desde el infinito».** Una fuente puntual en 2D es en
realidad una fuente lineal, así que un difusor o una barrera se interrogan
mejor con un frente plano. Dos herramientas lo cubren, cada una
unidireccional por su propio mecanismo:
- `sim.add_plane_wave(direction, center=..., width=..., wavelength=...)`
superpone un paquete gaussiano (opcionalmente con portadora sinusoidal)
como condición inicial viajando hacia `"down"`, `"up"`, `"left"` o
`"right"`; la velocidad consistente con el leapfrog escrita medio paso
atrás es lo que la hace unidireccional, y tras el frente la energía
residual queda al nivel del ruido numérico. Es lo que usa la animación
de difusión del QRD.
- `PlaneWaveSource(direction, waveform, offset=...)` registrada con
`add_source()` inyecta una onda plana sostenida sobre una línea de
celdas: la presión incidente y la velocidad de la cara adyacente se
excitan a la vez, de modo que el campo lanzado es plano en la
transversal a precisión de máquina y lo que se dispersa hacia atrás
cruza la línea intacto; con una esponja configurada detrás queda
absorbido.
```python
from phonometry.simulation import FDTD2D, CWSource, PlaneWaveSource
sim = FDTD2D(343.0, 0.01, shape=(160, 80), sponge_width=20,
sponge_sides=("top", "bottom"))
tone = CWSource(0, 0, frequency=1000.0) # reutilizada como forma de onda
sim.add_source(PlaneWaveSource("down", tone.value, offset=22))
# o, para un paquete único: sim.add_plane_wave("down", center=0.4,
# width=0.08, wavelength=0.34)
```
Las dos afirmaciones cuantitativas de ese punto son medibles, y la figura de
abajo las mide sobre exactamente esa escena. «Plano en la transversal a
precisión de máquina» es literal: pasado el transitorio de llenado, cada
columna del campo asentado lleva el *mismo* valor float64, así que la mayor
diferencia a lo ancho de una fila es 0.0 y no simplemente pequeña.
«Unidireccional» no lo es: la línea de inyección deja escapar un poco hacia
atrás, y lo que queda detrás solo es pequeño porque ahí está la esponja para
comérselo: 1,4 × 10⁻⁴ de la energía del campo, es decir −38,4 dB, se queda
en las 20 filas de esponja de detrás de la línea, pero justo detrás de la
línea la presión está apenas unos 26 dB por debajo de la onda que avanza.
Aleja la línea de su esponja, u olvídate de la esponja, y esos 26 dB son lo
que se te viene de vuelta.
*Qué aspecto tiene un lanzador unidireccional bien configurado: el frente es
plano hasta el último bit de un float64 (centro), el campo que avanza
mantiene amplitud unidad, y el residuo de detrás de la línea de inyección
arranca 26 dB por debajo y queda enterrado por la esponja que lleva ese lado
(derecha). La esponja no es un remate opcional: sin ella el residuo del
reverso se refleja en el borde superior y vuelve a atravesar la medición.*
Mostrar el código de esta figura
```python
# `FDTD2D`, `CWSource` y `PlaneWaveSource` son los nombres importados arriba.
sim = FDTD2D(343.0, 0.01, shape=(160, 80), sponge_width=20,
sponge_sides=("top", "bottom"))
sim.add_source(PlaneWaveSource(
"down", CWSource(0, 0, frequency=1000.0).value, offset=22))
sim.run(700) # rampa + llenado + unos periodos
body = sim.p[60:140, :]
print(float(np.abs(np.diff(body, axis=1)).max())) # 0.0 exactamente plano
back = float((sim.p[:20, :] ** 2).sum()) / float((sim.p ** 2).sum())
print(round(10 * np.log10(back), 1)) # -38.4 dB de energía
rms = np.sqrt((sim.p ** 2).mean(axis=1))
fig, (ax_f, ax_t, ax_l) = plt.subplots(1, 3, figsize=(13.5, 5.0))
ax_f.imshow(sim.p, cmap="RdBu_r", vmin=-1.05, vmax=1.05)
ax_t.plot(sim.p[80, :]) # corte transversal
ax_l.plot(20 * np.log10(rms / rms[60:140].mean()), np.arange(160))
plt.show()
```
La geometría se **rasteriza**: `obstacle_mask` marca celdas rígidas, y toda
cara en contacto con una celda marcada se cierra (de nuevo la Ec. 4.32), así
que paredes, barreras y dispersores de cualquier forma son simples arrays
booleanos. Una cara está abierta o cerrada, sin nada intermedio, así que una
superficie que corre a lo largo de un eje de la malla queda representada
exactamente y cualquier otra se convierte en una **escalera** cuyos peldaños
miden una celda. Esos peldaños dispersan energía que la superficie real no
dispersa, desplazan la superficie efectiva hasta media celda y, como el error
es una fracción fija de celda y no de longitud de onda, no se encoge al
alargarse la onda. Es peor a incidencia rasante y en lo alto de la banda
resuelta. La regla de trabajo es que un reflector inclinado o curvo necesita
una malla más fina de la que pediría por sí sola la regla de dispersión de
las diez celdas, con el mismo diagnóstico que para la dispersión: reduce
$\Delta x$ a la mitad y confirma que el campo dispersado converge en vez de
limitarse a cambiar. Por eso el panel mallado de la sección 4 corre a medio
milímetro, mucho más fino de lo que exige su longitud de onda de 17 cm: allí
la resolución la fija la geometría, no la longitud de onda.
Cada lado del dominio puede llevar su propia condición de contorno:
- ``"rigid"`` (por defecto): un reflector perfecto, $R = +1$.
- ``"absorbing"``: una capa esponja de `absorbing_layer_cells` celdas cuya
tasa de absorción crece cuadráticamente, emulando un contorno abierto (el
precursor sencillo de las capas perfectamente adaptadas de la sección
4.2.3 de Attenborough y Van Renterghem).
- una **impedancia específica real** $Z$ en Pa·s/m (un escalar o un valor
por celda de borde): el contorno de reacción local de las Ecs. 4.33-4.35,
actualizado implícitamente, con coeficiente de reflexión a incidencia
normal $R = (Z - \rho c)/(Z + \rho c)$; $Z = \rho c$ es anecoico.
**Las dos prestaciones no son intercambiables, y solo una de ellas es
blanda.** Los cuatro *lados del dominio* pueden ser rígidos, esponja o una
impedancia real de reacción local; todo lo dibujado en `obstacle_mask` es
perfectamente rígido por construcción. Una sala modelada como paredes
interiores dentro de un dominio mayor tiene, por tanto, paredes duras hagan
lo que hagan los bordes, y no hay forma de darle a una superficie interior un
coeficiente de absorción. Dos vías ablandan una superficie: ponerla en un
borde del dominio y darle a ese borde una impedancia, o respaldarla con una
región con pérdidas mediante el mapa `damping` de más abajo. Cuando basta con
un borde, la conversión desde un coeficiente de absorción objetivo a
incidencia normal son dos pasos: la magnitud de la reflexión es
$|R| = \sqrt{1-\alpha}$, y para una impedancia real
$Z = \rho c\,(1+R)/(1-R)$, así que $\alpha = 0{,}5$ pide
$Z \approx 5{,}83\,\rho c$ y $\alpha = 0{,}9$ alrededor de $1{,}92\,\rho c$. La
salvedad viene con ello: $R$ solo es exacto a incidencia normal, así que un
borde de reacción local absorbe menos cuanto más oblicua es la incidencia, y
una sala rodeada de bordes así decae más despacio de lo que sugiere su
$\alpha$ nominal.
**Pérdidas volumétricas: el mapa `damping`.** `damping` es una tasa de caída
de amplitud $\sigma$ por celda, en s⁻¹, aplicada a la vez a la presión y a
las dos componentes de la velocidad. Como ambos campos decaen juntos, una
onda plana dentro de una región uniforme con pérdidas sigue
$k = (\omega - j\sigma)/c$: la amplitud cae exponencialmente en el espacio a
$\sigma/c$ neperios por metro mientras la impedancia característica se
mantiene real en $\rho c$. Eso es un *fluido equivalente* (misma velocidad
del sonido, misma densidad, pérdidas independientes de la frecuencia) y es
la única pérdida volumétrica que tiene el esquema, de modo que también es la
única manera de hacer que absorba una superficie *interior*. Dos vías para
llegar a un valor. Como sustituto de la absorción de una sala,
$\sigma = 6{,}91/T_{60}$ produce una caída reverberante de $T_{60}$ segundos,
así que una caída de 0,5 s son $\sigma = 13{,}8$ s⁻¹. Como probeta, la pérdida
por metro son $\sigma/c$ neperios ($8{,}69\,\sigma/c$ dB/m), y el valor se
ajusta hasta que la muestra modelada reproduce una absorción medida. Los dos
límites merecen decirse sin rodeos: como la pérdida es independiente de la
frecuencia y la impedancia se mantiene real, esto **no** es un modelo de
Delany-Bazley ni de Johnson-Champoux-Allard y no reproducirá la dependencia
con la frecuencia de un absorbente real; y `FDTD2D` acepta un escalar *o* un
mapa `(ny, nx)`, mientras que la función de conveniencia `fdtd_simulation`
solo admite un escalar, así que un absorbente región a región hay que
construirlo sobre el motor.
```python
# Un conducto con pérdidas de tres filas: la caída por metro es todo el
# modelo. (`simulation` es el nombre que importó la sección 1.)
sigma = 200.0 # tasa de caída de amplitud [1/s]
lossy = simulation.FDTD2D(343.0, 0.01, shape=(3, 900), damping=sigma)
x = (np.arange(900) + 0.5) * 0.01
lossy.p[:] = np.exp(-(((x - 1.0) / 0.15) ** 2))[None, :]
peaks = []
for _ in range(round(0.020 / lossy.dt)):
lossy.step()
peaks.append((lossy.p[1, 200], lossy.p[1, 500]))
peaks = np.abs(np.asarray(peaks))
print(round(float(peaks[:, 1].max() / peaks[:, 0].max()), 3)) # 0.173
print(round(float(np.exp(-sigma * 3.0 / 343.0)), 3)) # 0.174 exacto
print(round(8.686 * sigma / 343.0, 2)) # 5.06 dB por metro
print(round(6.91 / 13.82, 2)) # 0.5 s de T60 con sigma = 13.82
```
La vía calibrada está al final de esta sección: la guía del tubo de
impedancia excita este mismo esquema con un mapa `damping` y pasa el
resultado por la reducción de ISO 10534-2, de modo que la absorción de la
muestra modelada sale de una medición virtual y no del valor que entró.
El motor de avance `FDTD2D` también es público: expone `step()`, `run()`,
los arrays del campo y la energía, para quien necesite acceso fotograma a
fotograma (las animaciones de la documentación lo usan directamente). Un
pulso plano lanzado por un conducto contra un borde de impedancia reproduce
el coeficiente de reflexión de libro:
```python
from phonometry import simulation
rho, c, dx = 1.2, 343.0, 0.01
sim = simulation.FDTD2D(c, dx, rho=rho, shape=(3, 1200),
edge_impedance={"right": 3.0 * rho * c})
x = (np.arange(1200) + 0.5) * dx
sim.p[:] = np.exp(-(((x - 6.0) / 0.15) ** 2))[None, :] # pulso plano
trace = []
for _ in range(int(round(0.032 / sim.dt))):
sim.step()
trace.append(sim.p[1, 900])
trace = np.asarray(trace)
t = (np.arange(trace.size) + 1) * sim.dt
t_return = 6.0 / c + 3.0 / c # vía la pared, de vuelta a x = 9 m
incident = trace[t < t_return - 0.001].max()
echo = trace[t > t_return]
print(round(float(echo[np.abs(echo).argmax()] / incident), 2)) # 0.5
# (Z - rho c)/(Z + rho c) = (3 - 1)/(3 + 1) = +0.5
```
Antes de gastar un solo paso de tiempo, `sim.plot_geometry()` dibuja el
dominio configurado: bordes, esponjas, obstáculos, fuentes y las sondas que
piensas registrar. Pillar una esponja en el lado equivocado o una sonda mal
colocada cuesta segundos aquí y una ejecución entera después.
*Todo lo que la simulación verá, antes de correrla: las capas de esponja se
comen los contornos izquierdo y derecho, el borde superior lleva la
impedancia anecoica $\rho c$, el borde inferior sin tratar queda rígido, y las
dos sondas quedan libres del obstáculo y de las esponjas.*
Mostrar el código de esta figura
```python
from phonometry import simulation
mask = np.zeros((60, 90), dtype=bool)
mask[25:35, 40:44] = True
sim = simulation.FDTD2D(343.0, 0.05, shape=(60, 90), sponge_width=8,
sponge_sides=("left", "right"),
edge_impedance={"top": 413.0}, obstacle_mask=mask)
sim.add_source(simulation.GaussianPulse(10, 30, width=1e-3))
# Comprueba el dominio antes de correrlo: aún no se ha dado ningún paso.
sim.plot_geometry(probes=[(3.0, 1.5), (4.0, 2.0)], language="es")
plt.show()
```
### Comprobaciones de que una ejecución es utilizable
`plot_geometry` pilla una escena mal construida. Estas seis pillan una
ejecución mal *calculada*, y cada una de ellas es una línea o dos:
1. **La energía no debe crecer con el tiempo.** `FDTD2D.energy()` devuelve la
energía total del campo en julios por metro de profundidad. Sube mientras
la fuente inyecta y luego se asienta en una meseta, con una fluctuación de
en torno al 1 % porque la presión y la velocidad están separadas medio
paso temporal; en un dominio cerrado sin pérdidas la meseta es plana, y con
una esponja o un mapa `damping` decae. Un total que sube después de apagar
la fuente es energía que se está creando, y la ejecución no sirve por
verosímil que parezca la instantánea.
2. **`np.isfinite(sim.p).all()`** pilla una divergencia antes de que envenene
una FFT en silencio: en cuanto una sola celda es `inf`, todo espectro
sacado de ese registro es `nan` y el gráfico sale vacío en vez de salir
mal.
3. **Ninguna sonda dentro de una esponja, de un obstáculo o de una celda de
fuente.** Las tres leen algo que no es el campo: una celda de obstáculo
vale exactamente 0.0 durante toda la ejecución, una celda de esponja lee
el campo ya atenuado camino de la salida, y una celda de fuente blanda lee
sobre todo la inyección; en la escena de la barrera de la sección 5, una
sonda sobre la celda de la fuente marca un pico 6,7 veces mayor que la
sonda A, a 0,4 m, y una sonda dentro de la barrera lee exactamente 0.0
durante toda la ejecución.
4. **Aplica la regla de las diez celdas con la velocidad del sonido *más
baja* del mapa**, a la frecuencia más alta que el análisis vaya a usar de
verdad, no con el valor del aire ni con la frecuencia más alta que la
fuente excite por casualidad (la sección 3 dimensiona ambas).
5. **Vuelve a ejecutar con la mitad de tamaño de celda y confirma que la
respuesta se mueve menos que la tolerancia que necesita tu afirmación.** El
esquema es de segundo orden, así que una magnitud limitada por la
dispersión debería moverse unas cuatro veces menos al partir por la mitad
el paso de malla; una magnitud que se mueve otro tanto está limitada por
otra cosa: el escalonado de arriba, o la longitud del registro.
6. **En una escena que se pretende anecoica, ejecútala una vez sin el
dispersor** y comprueba que el residuo queda por debajo del nivel que
piensas dar. Esa resta es exactamente lo que automatiza
`ContourPhasors.subtract()` en la sección 4, y lo que hace la medición de
absorción in situ de ISO 13472-1 sobre una carretera real.
Un dominio de tres filas con paredes rígidas es un tubo de onda plana, y con
el mapa de `damping` por celda una muestra porosa se convierte en un fluido
equivalente, lo que convierte a la propia simulación en una probeta medible: la
[guía del tubo de
impedancia](/phonometry/es/materials/absorbers/impedance-tube/) ejecuta
las mediciones ISO 10534-2 y ASTM E2611 de forma virtual sobre exactamente
ese dominio, animaciones incluidas, y recupera la absorción y la pérdida por
transmisión analíticas de la muestra modelada con las propias cadenas de
reducción de la biblioteca (los contrastes de `tests/simulation` lo ejecutan
en cada commit).
## 3. Dimensionar una ejecución
La sección 6 deduce la única regla que da la mayoría de los textos de FDTD
(diez celdas por longitud de onda más corta) y esa regla decide solo
$\Delta x$. Otros cuatro números deciden si la ejecución dice algo en
absoluto: el **ancho de banda de la fuente**, la **duración de la
ejecución**, el **grosor de la esponja** y, para una medición en régimen
estacionario, el **transitorio** que se tira antes de la **ventana DFT** que
se conserva. Cada uno tiene su regla, y cada uno falla en silencio cuando
está mal.
**Ancho de banda de la fuente.** Un `GaussianPulse` de semianchura `width`
es $s(t) = A\,e^{-((t-t_0)/\text{width})^2}$, cuyo espectro de magnitud es
$e^{-(\pi f\,\text{width})^2}$: queda 20 dB por debajo en
$f = 0{,}483/\text{width}$ y 40 dB por debajo en $0{,}683/\text{width}$.
Trabaja en la dirección del lector: elige $f_\text{max}$, saca de ahí
$\Delta x$, y toma después $\text{width} \ge 0{,}5/f_\text{max}$; así el
pulso está al menos 20 dB por debajo donde la malla deja de resolver. De
forma equivalente, en términos de la malla,
$\text{width} \ge 4{,}83\,\Delta x/c_\text{min}$, que con $\Delta x = 1$ cm
en aire son $1{,}4 \times 10^{-4}$ s, y es la razón por la que los ejemplos
de aquí usan $3 \times 10^{-4}$ s, excitando hasta unos 1,6 kHz dentro de
los 3,4 kHz que la malla resuelve. La energía por encima de $f_\text{max}$
no es un error en sí misma; simplemente viaja a la velocidad equivocada y
vuelve como una cola de rizado detrás del pulso. La ejecución en caja rígida
de la sección 6 lee un modo a 299 Hz sobre una malla cuyo límite de diez
celdas está en 1,7 kHz, así que su pulso más estrecho de
$2 \times 10^{-4}$ s no le cuesta nada. Para un estudio a una sola
frecuencia, `CWSource` esquiva la cuestión por completo.
**Duración.** La fijan tres cosas, y manda la mayor. *Geometría*: como
mínimo el tiempo de cruce del dominio $\sqrt{L_x^2 + L_y^2}/c$, más un viaje
de ida y vuelta al reflector por cada eco que importe. *Decaimiento*: lo
bastante larga para que el rasgo se vea por encima de lo que aún resuena.
*Espectro*: una FFT de sonda resuelve $1/T$, así que separar dos modos que
distan 3 Hz cuesta un tercio de segundo de registro por fina que sea la
malla, que es exactamente por lo que la ejecución de modos de sala de abajo
dura 0,35 s y la de la barrera 9 ms.
**Grosor de la esponja.** La capa hace crecer cuadráticamente una tasa de
absorción $\sigma$ desde cero en su cara interior hasta un máximo en el
borde exterior. Como $\sigma$ hace decaer a la vez la presión y la
velocidad, la impedancia característica dentro de ella se mantiene en
$\rho c$ y una región con pérdidas *uniforme* no reflejaría nada; todo lo
que vuelve lo refleja el *gradiente*. Eso convierte la regla en una regla en
**celdas**, no en longitudes de onda, y merece la pena ver la medición,
tomada a incidencia normal sobre un conducto de onda plana:
| `sponge_width` | 5 | 10 | 20 | 40 | 60 |
|---|---|---|---|---|---|
| eco/incidente a 1 kHz | −36 dB | −49 dB | −62 dB | −76 dB | −85 dB |
| eco/incidente a 125 Hz | — | — | −61 dB | — | — |
Una capa de 20 celdas devuelve −61 dB a 125 Hz, donde mide 0,07 longitudes
de onda de grosor, y −64 dB a 2 kHz, donde mide 1,2 longitudes de onda: el
residuo es plano dentro de 3 dB sobre un margen de frecuencias de 16:1.
Veinte celdas son un buen suelo; cuarenta, generosas. `sponge_reflection`
(por defecto $10^{-4}$) fija la amplitud de ida y vuelta para la que la
rampa está *diseñada*, y no es una promesa: con 20 celdas el residuo medido
es de −51 dB para $10^{-2}$, −62 dB para $10^{-4}$ y −59 dB para $10^{-6}$,
porque un objetivo más agresivo hace la rampa más abrupta y una rampa más
abrupta refleja más. El valor por defecto está cerca del óptimo, y la manera
de comprar más es el grosor.
Dos cosas que la esponja no hace. La primera, es mucho más débil a
incidencia **oblicua y rasante**, donde una onda que corre a lo largo de la
capa atraviesa muy pocas celdas de la rampa por longitud de onda recorrida.
Toma esa misma capa de 20 celdas, pon una fuente de pulso 1 m por encima y
mide lo que la capa añade a una sonda a la misma altura volviendo a ejecutar
con la capa a 5 m y restando: lo que añade es −30 dB respecto a la llegada
directa a 27° de la normal, −18 dB a 45°, −10 dB a 60° y −3 dB a 80°, frente
a los −62 dB que esa misma capa da a incidencia normal. Un dominio largo y
poco profundo forrado de esponjas es, por tanto, un conducto y no un campo
libre. La prueba honesta es la que se acaba de describir: aleja el contorno
y confirma que la traza de la sonda no cambia. La segunda, la esponja no es
sitio para colocarse: sondas, fuentes y obstáculos van todos fuera de la
capa (comprobación 3 de arriba).
**Transitorio y ventana.** Una medición en régimen estacionario se ejecuta
en dos partes: llenar el dominio, llamar a `probe.reset()` e integrar
después sobre un número exacto de periodos de la fuente. El transitorio que
se descarta es el tiempo de llenado del dominio $\sqrt{L_x^2+L_y^2}/c$, más
una travesía por cada reflexión que todavía importe, más el arranque de la
fuente (`ramp_cycles / frequency`, tres periodos por defecto), más el tiempo
que tarde en excitarse cualquier geometría resonante, del orden de $Q$
periodos, que es toda la razón por la que la ejecución del metadifusor de
abajo espera 8 ms donde el monopolo desnudo espera 4,5 ms: sus cuellos y
cavidades son resonadores y tienen que llenarse. La ventana que viene
después debería abarcar un número entero de periodos de la fuente, para que
la DFT al vuelo no vea truncamiento; con diez periodos sobra.
En una línea: $\Delta x$ a partir de $f_\text{max}$ y $c_\text{min}$,
$\Delta t$ a partir del número de Courant y $c_\text{max}$, `width` a partir
de $\Delta x$, la esponja a partir del presupuesto de celdas, la duración a
partir de la geometría y de la resolución espectral, el transitorio a partir
del dominio y de sus resonancias, y la ventana a partir del periodo de la
fuente.
## 4. Del campo cercano al campo lejano
Una respuesta polar o un coeficiente de difusión son magnitudes de campo
lejano, pero una caja FDTD termina a un par de longitudes de onda del
dispersor. La pareja `add_contour_probe` / `far_field_from_contour` salva
esa distancia con la integral de Kirchhoff-Helmholtz 2D: la sonda pliega la
presión y la velocidad normal saliente del régimen estacionario, sobre un
rectángulo cerrado de caras de celda, en acumuladores complejos (una DFT al
vuelo por punto y frecuencia, así que un régimen de onda continua no
almacena historias temporales), y la integral propaga esos fasores hasta el
infinito con la función de Green de espacio libre $-(j/4)\,H_0^{(2)}(kR)$, la
misma construcción con la que los códigos FEM de onda completa obtienen los
patrones de dispersión a partir del campo cercano. Dos propiedades
sostienen el esquema: la malla escalonada ya guarda $v_n$ exactamente en
las caras del contorno, y todo campo cuyas fuentes queden *fuera* del
contorno integra a nada (extinción), de modo que los fasores del campo
total de una dispersión con onda plana se transforman directamente en el
campo lejano *dispersado*. En el esquema discreto esa cancelación deja
un residuo de dispersión numérica de malla (por debajo de 0,01 dB en las
escenas validadas aquí); `ContourPhasors.subtract()` lo elimina con una
pasada de referencia sin dispersor cuando esa última fracción importa.
Nada de eso se ve en la aritmética de índices que hacen los fragmentos, así
que aquí está la escena que construye de verdad la ejecución del metadifusor
de abajo, en su proporción real.
```python
from phonometry import simulation
c0, dx, f0 = 343.0, 0.005, 2000.0
sim = simulation.FDTD2D(c0, dx, shape=(300, 300), sponge_width=40)
sim.add_source(simulation.CWSource(ix=150, iy=150, frequency=f0))
probe = sim.add_contour_probe(90, 210, 90, 210, frequencies=[f0])
sim.run(round(4.5e-3 / sim.dt)) # agotar el transitorio
probe.reset() # e integrar la ventana DFT
sim.run(round(10.0 / f0 / sim.dt))
pattern = simulation.far_field_from_contour(
probe.phasors(f0), np.arange(0.0, 360.0, 5.0),
origin=(150.5 * dx, 150.5 * dx))
levels = 20 * np.log10(np.abs(pattern))
print(round(float(levels.max() - levels.min()), 2)) # 0.04 dB de rizado:
# una fuente lineal es omnidireccional, y su nivel coincide con la
# función de Green 2D de espacio libre en 0,11 dB
```
Los oráculos analíticos tras esos números se ejecutan en cada commit
(`tests/simulation`): el patrón del monopolo reconstruido es plano dentro
de 0,05 dB y se asienta sobre el nivel de la función de Green 2D dentro de
0,11 dB, una pareja en contrafase reproduce el factor de array de dos
fuentes dentro del 0,4 % del pico, y un contorno que no encierra la fuente
se transforma en menos del 1 % de uno que sí. Sobre ellos se apoyan dos
contrastes de panel mallado a 2 kHz: un difusor de residuo cuadrático con
pozos de hasta 27,4 cm, cuya respuesta polar NTFF sigue la predicción de
Fraunhofer de `predict_diffuser_polar_response` (correlación de patrones
de 0,94, coeficiente de difusión direccional ISO 17497-2 dentro de
0,08), y el metadifusor en sublongitud de onda profunda de abajo.
La cadena de campo lejano es la que por fin cierra de extremo a extremo el
círculo del metadifusor: el panel de la Tabla 1 del artículo, con cada
ranura, cuello y cavidad mallados a 0,5 mm, llevado al régimen
estacionario por una onda plana, capturado sobre un contorno y
transformado, frente al modelo
[`metadiffuser_polar_response`](/phonometry/es/reference/api/materials/metadiffuser/)
(TMM + Fraunhofer) del módulo de materiales. Es la réplica en la biblioteca
del contraste TMM frente a FEM que documenta el propio artículo de los
metadifusores, pequeñas discrepancias incluidas: el modelo de matrices de
transferencia homogeneiza cada celda de 7 cm en un único coeficiente de
reflexión de reacción local e ignora el acoplamiento evanescente entre
bocas vecinas, así que la estructura de lóbulos coincide mientras algunos
nulos se desplazan unos grados.
*Dos rutas independientes al mismo campo lejano: el panel mallado por
onda completa en el dominio del tiempo (continuo) frente a
la cadena de matrices de transferencia homogeneizada (discontinuo). Los
lóbulos coinciden; los nulos, sensibles a las correcciones de extremo
milimétricas que el TMM modela analíticamente, se desplazan unos grados.*
Mostrar el código de esta figura
```python
from phonometry import (HelmholtzResonator, MetadiffuserWell,
metadiffuser_polar_response, simulation)
c0, dx, f0, pitch = 343.0, 0.0005, 2000.0, 0.07
rows = [(14.7, 13.0, 16.4, 6.2, 9.0), (30.9, 9.1, 4.3, 3.5, 9.0),
(30.9, 9.1, 4.3, 3.5, 9.0), (15.7, 13.3, 17.0, 6.3, 9.0),
(20.3, 18.0, 20.7, 3.2, 9.0)] # Tabla 1: h, l_n, l_c, w_n, w_c
# Mallar el panel real: una losa rígida con las ranuras, cuellos y
# cavidades excavados a 0,5 mm (el cuello más estrecho, 3,2 mm, abarca
# seis celdas).
sponge, gap, marg, front = 60, 60, 20, 40
face = round(5 * pitch / dx)
lat = marg + gap + sponge
r_face = sponge + gap + front
slab = round(0.023 / dx) # panel de 2 cm + fondo de 3 mm
mask = np.zeros((r_face + slab + marg + gap + sponge,
face + 2 * lat), dtype=bool)
mask[r_face:r_face + slab, lat:lat + face] = True
for n, (h, ln, lc, wn, wc) in enumerate(rows):
xs = (n + 0.12) * pitch
c0s, c1s = lat + round(xs / dx), lat + round((xs + h * 1e-3) / dx)
mask[r_face:r_face + round(0.02 / dx), c0s:c1s] = False
for m in range(2): # dos resonadores por ranura
ym, xn = (m + 0.5) * 0.01, xs + h * 1e-3
r0 = r_face + round((ym - 0.5e-3 * wn) / dx)
r1 = r_face + round((ym + 0.5e-3 * wn) / dx)
mask[r0:r1, c1s:lat + round((xn + ln * 1e-3) / dx)] = False
r0 = r_face + round((ym - 0.5e-3 * wc) / dx)
r1 = r_face + round((ym + 0.5e-3 * wc) / dx)
mask[r0:r1, lat + round((xn + ln * 1e-3) / dx):
lat + round((xn + (ln + lc) * 1e-3) / dx)] = False
sim = simulation.FDTD2D(c0, dx, shape=mask.shape, sponge_width=sponge,
cfl=0.9, obstacle_mask=mask) # 340 celdas/lambda:
sim.add_source(simulation.PlaneWaveSource( # dispersión despreciable
"down", simulation.CWSource(0, 0, f0).value, offset=sponge))
probe = sim.add_contour_probe(lat - marg, lat + face + marg - 1,
r_face - front, r_face + slab + marg - 1,
frequencies=[f0])
sim.run(round(8e-3 / sim.dt)) # agotar transitorio y resonadores
probe.reset()
sim.run(round(10.0 / f0 / sim.dt)) # ventana DFT de 10 periodos (en pasos)
angles = np.arange(-90.0, 90.1, 5.0) # desde la normal del panel
pattern = simulation.far_field_from_contour(
probe.phasors(f0), angles - 90.0, # la normal apunta según -y
origin=((lat + face / 2.0) * dx, r_face * dx))
levels = 20 * np.log10(np.abs(pattern) / np.abs(pattern).max())
wells = [MetadiffuserWell(h * 1e-3,
(HelmholtzResonator(ln * 1e-3, wn * 1e-3,
lc * 1e-3, wc * 1e-3),) * 2)
for h, ln, lc, wn, wc in rows]
model = metadiffuser_polar_response(f0, wells, depth=0.02, period=pitch,
angles=angles, periods=1)
ax = model.plot(color="#1f77b4", marker="", linestyle="--",
label="Modelo TMM + Fraunhofer", language="es")
ax.plot(np.radians(angles), levels, color="#d62728", lw=2.2,
label="FDTD + NTFF, panel mallado a 0,5 mm")
ax.set_ylim(-40.0, 2.0)
ax.legend(loc="lower center")
plt.show()
```
Por qué se desplazan los nulos es una pregunta sobre geometría, y merece la
pena poner las dos geometrías una al lado de la otra: el panel que el modelo
tasa y el panel sobre el que avanza el esquema.
*Lo que homogeneiza la matriz de transferencia, frente a lo que ve de verdad
el esquema. El modelo convierte cada celda de 70 mm en un único coeficiente
de reflexión de reacción local; la máscara conserva cada ranura, cuello y
cavidad, y el acoplamiento evanescente entre bocas vecinas sale gratis
porque las bocas están realmente ahí. El cuello más estrecho mide 3,2 mm
(seis celdas con $\Delta x = 0{,}5$ mm) y es él, y no la longitud de onda de
17,2 cm, quien fuerza la malla. En esa diferencia está el origen de los
nulos desplazados.*
Mostrar el código de esta figura
```python
# La vista del modelo (panel superior) la dibuja la propia biblioteca, a
# partir de los mismos `wells` que construyó el bloque de arriba:
# from phonometry.materials import plot_metadiffuser_panel_geometry
# plot_metadiffuser_panel_geometry(wells, depth=0.02, period=0.07)
# La vista del esquema: la máscara booleana que malló el bloque de arriba, y
# una celda suya ampliada.
fig, (ax_all, ax_cell) = plt.subplots(2, 1, figsize=(11, 6))
ax_all.imshow(mask[r_face - 40:r_face + slab + 20, lat - 20:lat + face + 20],
cmap="Greys", interpolation="nearest")
cell = lat + round(4 * pitch / dx)
ax_cell.imshow(mask[r_face - 40:r_face + slab + 20,
cell:cell + round(pitch / dx)],
cmap="Greys", interpolation="nearest")
plt.show()
```
## 5. Cuándo merece la pena una simulación de ondas, y qué cambia el 2D
FDTD amortiza su coste cuando la **geometría gobierna la física**: difracción
en torno a una barrera o a través de una abertura, interferencia de caminos
directo y reflejado, dispersión por obstáculos, comportamiento modal de
recintos de forma irregular, refracción en un gradiente de velocidad del
sonido. Una sola ejecución captura **todas las frecuencias a la vez** (un
pulso excita toda la banda; una FFT de una sonda da el espectro), mientras
que un método en el dominio de la frecuencia necesita resolver el sistema
una vez por cada frecuencia.
La ejecución de la barrera de abajo es el arquetipo. Nada en su entrada dice
«difracción»: la barrera son cuatro columnas de un array booleano, la fuente
es un pulso, y la llegada a la zona de sombra está ahí de todas formas.
*Una ejecución, dos vistas. La instantánea (izquierda, a los 6,5 ms) recoge
el frente directo, la reflexión que vuelve pasando junto a la fuente y la
onda que dobla el borde superior de la barrera. En las trazas (derecha), la
sonda A, con visión directa, registra el pulso directo y después la
reflexión de la barrera, mientras que la sonda B, en sombra, solo registra
la llegada más débil y más tardía que vino rodeando el borde: una amplitud y
un retardo que ninguna fórmula de divergencia más atenuación produce por sí
sola.*
Mostrar el código de esta figura
```python
from phonometry import simulation
# Un campo libre de 3.0 x 2.0 m (bordes absorbentes) con una barrera
# rígida delgada: la sonda A ve el pulso directo más la reflexión de la
# barrera, la sonda B queda en sombra y solo recibe la onda difractada.
mask = np.zeros((200, 300), dtype=bool)
mask[60:, 150:154] = True
res = simulation.fdtd_simulation(
343.0, 0.01, 9.0e-3, shape=(200, 300),
sources=[simulation.GaussianPulse(ix=60, iy=100, width=3.0e-4)],
probes=[(100, 100), (240, 100)],
obstacle_mask=mask,
boundaries="absorbing", absorbing_layer_cells=30,
snapshot_every=75,
)
fig, (ax_f, ax_p) = plt.subplots(
1, 2, figsize=(12.5, 5.0), gridspec_kw={"width_ratios": [1.25, 1.0]})
res.plot(kind="snapshot", frame=7, ax=ax_f, language="es")
res.plot(ax=ax_p, language="es")
plt.tight_layout()
plt.show()
```
A la inversa, cuando existe una forma cerrada validada (reverberación
estadística, atenuación en exteriores de ISO 9613-2, fuentes imagen en una
sala rectangular), esa forma cerrada es miles de veces más barata, y esta
simulación es el contraste y el demostrador más que el sustituto. La sección
6 muestra uno de esos contrastes reproduciendo los modos analíticos de la
sala.
**Lo que cuesta una ejecución.** El trabajo son celdas × pasos, y el avance
monohilo en float64 hace del orden de $10^8$ actualizaciones de celda por
segundo (medido: $8 \times 10^7$ en la máquina que renderiza estas figuras).
La escena de la barrera de arriba son 200 × 300 celdas durante 728 pasos:
$4{,}4 \times 10^7$ actualizaciones, medio segundo. La escena del
metadifusor de la sección 4 se malla a 0,5 mm en 346 × 980 celdas y corre
unos 14 000 pasos: $4{,}8 \times 10^9$ actualizaciones, alrededor de un
minuto. El escalado entre esas dos es la parte que conviene recordar. En 2D,
reducir $\Delta x$ a la mitad cuadruplica el número de celdas *y*, por la
condición de Courant, duplica el número de pasos, así que el trabajo se
multiplica **por ocho** a cambio de dividir entre cuatro el error de
dispersión. La regla de las diez celdas es, por tanto, un suelo en el que
asentarse y no un objetivo que batir.
La memoria sigue la misma aritmética. El esquema guarda doce mapas
`float64` (la presión, las dos componentes de la velocidad, la velocidad del
sonido, la densidad, el módulo de compresibilidad, las densidades en las
caras y los factores de caída): unos 96 bytes por celda, 113 con una máscara
de obstáculos, así que la escena del metadifusor vive en unos 39 MB. Las
instantáneas almacenadas son el mando que nadie nombra: cada una son otros 8
bytes por celda, así que `snapshot_every=1` sobre esa escena pediría del
orden de 38 GB. Elige la cadencia a partir del número de fotogramas que
quieres de verdad; las animaciones de esta documentación usan unos
centenares. Las historias de sonda cuestan un float por sonda y por paso,
que no es nada. El orden práctico de trabajo es dimensionar la malla a
partir de $f_\text{max}$, estimar el número de pasos a partir de la
duración, y probar el conjunto una vez con el doble de tamaño de celda para
comprobar la geometría y la colocación de las sondas antes de pagar la
ejecución de verdad.
El dominio es **bidimensional**, y eso cambia la física, no solo el coste.
Una fuente puntual 2D es físicamente una **fuente lineal** infinita: su
amplitud se expande cilíndricamente como $1/\sqrt{r}$ (3,0 dB por
duplicación de distancia) en lugar del $1/r$ esférico (6,0 dB) de una
fuente puntual 3D, y la respuesta impulsional 2D arrastra una estela tras el
frente en vez de pasar limpiamente. Los *patrones* de interferencia y
difracción son fieles; los niveles absolutos y las tasas de caída no son los
de una sala 3D. Trata las ejecuciones 2D como secciones y demostraciones, y
valida cualquier afirmación cuantitativa 3D contra una forma cerrada o un
cálculo 3D.
## 6. Dispersión numérica y precisión
La malla discreta propaga cada frecuencia a una velocidad ligeramente
errónea: las longitudes de onda cortas se retrasan, así que un pulso agudo
desarrolla una cola de rizado y las resonancias se desplazan un poco. Esta
**dispersión numérica** es la contrapartida discreta de la Ec. 4.15; sobre
los ejes de una malla cuadrada la relación de dispersión del esquema es
$$
\sin\!\left(\frac{\omega\,\Delta t}{2}\right)
= \frac{c\,\Delta t}{\Delta x}\,
\sin\!\left(\frac{k\,\Delta x}{2}\right),
$$
con un error relativo de frecuencia a primer orden de magnitud
$(1 - S^2)\,(k\,\Delta x)^2 / 24$ sobre los ejes de la malla (la frecuencia
modelada queda por debajo de la real, así que el error con signo es
negativo), donde $S = c\,\Delta t / \Delta x$ es el número de Courant **por
eje**, que con celdas cuadradas es el valor de `cfl` dividido por $\sqrt{2}$:
el `cfl = 0.6` por defecto da $S = 0{,}424$ y $1 - S^2 = 0{,}82$; el error es
máximo exactamente sobre un
eje y se anula a lo largo de la diagonal de las celdas en el límite de
Courant $\mathrm{CN} = 1$. La regla práctica es resolver **al menos 10 celdas por
longitud de onda más corta**, $\Delta x \le c_\text{min} / (10 f_\text{max})$ con la velocidad
de sonido más baja del dominio: con exactamente 10 celdas la cota
$(k\,\Delta x)^2 / 24$ para números de Courant pequeños da alrededor del
1,6 %, que el factor $1 - S^2$ reduce a en torno al 1,4 % con el
``cfl = 0.6`` por defecto
(en un dominio heterogéneo el paso temporal lo fijan las celdas más
rápidas, así que las regiones lentas trabajan a un número de Courant local
menor y quedan más cerca de la cota del 1,6 %), y toda componente mejor
resuelta o fuera de eje es más precisa. Con
$\Delta x = 1$ cm el punto de 10 celdas en aire queda hacia los 3,4 kHz, y
reducir $\Delta x$ a la mitad divide el error entre cuatro (el esquema es de
segundo orden, y la batería de validación mide ese orden observado bajo
refinamiento de malla).
**Un pulso paga el triple.** El error de arriba es un error de velocidad de
*fase*: dice que un tono estacionario de número de onda $k$ viaja un poco
lento. Un pulso viaja a la velocidad de **grupo**, y derivar esa misma
relación de dispersión da
$v_\mathrm{g}/c = \cos\theta / \sqrt{1 - S^2\sin^2\theta}$ con
$\theta = k\,\Delta x / 2$, cuyo desarrollo es
$1 - (1 - S^2)(k\,\Delta x)^2/8$: la misma ley con un 8 en lugar de un 24,
así que un paquete de ondas llega tres veces más tarde de lo que sugiere la
regla de la fase. Con 10 celdas por longitud de onda y el `cfl` por defecto
eso es el 4,1 %, no el 1,4 %. Y como es un error de velocidad, se paga **por
metro recorrido**: la misma malla que va un 4 % lenta cuesta un centímetro
en un metro y un metro en veinticinco.
El clip de abajo mide las dos afirmaciones. Se lanza una sola ráfaga tonal
de 500 Hz por tres tubos planos idénticos en todo salvo en $\Delta x$
(137,2, 68,6 y 34,3 mm, es decir, 5, 10 y 20 celdas por longitud de onda), y
la respuesta continua exacta viaja en gris detrás de cada traza, así que el
retraso del paquete numérico es un hueco en pantalla y no un número en una
tabla.
*Aquí se leen tres cosas. Primera, el retraso crece: es de 0,14 m a los
2,4 ms y de 1,35 m al final, porque un error de velocidad se acumula con la
distancia. Segunda, al tubo de malla más gruesa le crece una **cola de
rizado por detrás** del paquete, no una reflexión por delante: las
componentes cortas de dentro de la ráfaga son las más lentas, así que llegan
las últimas. Tercera, los déficits medidos son del 17,2 %, 4,3 % y 1,1 %
frente al 17,3 %, 4,3 % y 1,1 % de la forma cerrada evaluada sobre el propio
espectro de la ráfaga, y los paquetes cruzan la línea de meta 3,3, 0,7 y
0,2 ms después que la onda exacta. Reducir $\Delta x$ a la mitad divide el
error entre cuatro, en las dos leyes y en la medición.*
La verificación funciona en los dos sentidos: aquí la forma cerrada ancla el
esquema, y una vez anclado es el esquema el que pone a prueba una forma
cerrada fuera de sus propias hipótesis. Una ejecución en caja rígida
reproduce los modos analíticos de la sala.
```python
from phonometry import simulation
lx, ly, dx = 1.0, 0.7, 0.02
nx, ny = round(lx / dx), round(ly / dx)
res = simulation.fdtd_simulation(
343.0, dx, 0.35, shape=(ny, nx),
sources=[simulation.GaussianPulse(ix=7, iy=5, width=2.0e-4)],
probes=[(nx - 4, ny - 3)],
)
p = res.pressures[0]
# Enventana el registro (termina con los modos todavía resonando sin
# amortiguar, y un truncamiento abrupto emborrona todos los picos) y
# rellénalo con ceros por ocho, lo que interpola la posición del pico sin
# añadir resolución real.
spec = np.abs(np.fft.rfft(p * np.hanning(p.size), n=8 * p.size))
freqs = np.fft.rfftfreq(8 * p.size, res.dt)
sel = (freqs > 250) & (freqs < 350)
print(round(0.5 * 343.0 * float(np.hypot(1 / lx, 1 / ly)), 1)) # 299.1 modo (1,1) exacto
print(round(float(freqs[sel][np.argmax(spec[sel])]), 1)) # 298.9 medido
```
*El espectro en la sonda del ensayo de caja rígida presenta sus picos sobre
las frecuencias modales analíticas
$f = (c/2)\sqrt{(n_x/L_x)^2 + (n_y/L_y)^2}$ (Kuttruff 6.ª ed., cap. 3). El
pico (1,1) cae un 0,05 % por debajo y la cota de arriba predice un 0,04 %,
pero esa coincidencia hay que leerla con cuidado: con 0,35 s la separación
bruta de bins es de 2,9 Hz y el relleno de ceros por ocho interpola hasta
0,36 Hz, así que el desplazamiento de 0,15 Hz está en el límite de lo que el
registro puede resolver. Establecer bien el error de dispersión pasa por
refinar la malla, no por mirar un pico con más ganas.*
Mostrar el código de esta figura
```python
from phonometry import simulation
lx, ly, dx, c = 1.0, 0.7, 0.02, 343.0
nx, ny = round(lx / dx), round(ly / dx)
res = simulation.fdtd_simulation(
c, dx, 0.35, shape=(ny, nx),
sources=[simulation.GaussianPulse(ix=7, iy=5, width=2.0e-4)],
probes=[(nx - 4, ny - 3)],
)
# Una línea: la historia de presión en la sonda de la que sale el espectro.
res.plot(language="es")
plt.show()
# La comprobación modal: espectro en la sonda frente a los modos analíticos.
p = res.pressures[0]
spec = np.abs(np.fft.rfft(p * np.hanning(p.size), n=8 * p.size))
freqs = np.fft.rfftfreq(8 * p.size, res.dt)
sel = (freqs >= 100) & (freqs <= 450)
fig, ax = plt.subplots()
ax.plot(freqs[sel], 20 * np.log10(spec[sel] / spec[sel].max()))
for mx, my in [(1, 0), (0, 1), (1, 1), (2, 0), (2, 1)]:
ax.axvline(0.5 * c * np.hypot(mx / lx, my / ly), ls=":", color="tab:red")
ax.set(xlabel="Frecuencia [Hz]", ylabel="Espectro en la sonda [dB re máx]",
ylim=(-60, 6))
plt.show()
```
Los tests anclan el esquema a oráculos analíticos de la misma manera:
autofrecuencias de caja y de conducto, tiempos de llegada en campo libre y
caída cilíndrica, el eco de la fuente imagen de una pared rígida, el
coeficiente de reflexión de impedancia de la sección 2 y la propia relación
de dispersión.
El `FDTDResult` congelado lleva el eje temporal, las historias de presión por
sonda, las posiciones de las sondas en metros, los metadatos de la malla, las
fuentes, las instantáneas opcionales del campo con sus tiempos y la máscara
de obstáculos; su `.plot()` dibuja las historias de las sondas, y
`.plot(kind="snapshot")` representa un campo registrado con la geometría
superpuesta.
La misma malla escalonada se extiende más allá de los fluidos: el
compañero `elastic_fdtd_simulation` integra sobre ella el sistema
velocidad-esfuerzo P-SV de Virieux (1986), añadiendo ondas de cizalla,
superficies libres por imagen de esfuerzos con ondas de Rayleigh, y el
acoplamiento fluido-sólido con conversión de modo, ondas de interfase de
Scholte y transmisión de placas sumergidas. Ese esquema tiene su
propia guía,
[Ondas elásticas y acoplamiento fluido-sólido](/phonometry/es/simulation/elastic-waves/).
## Qué cubre esta guía
El esquema FDTD presión-velocidad escalonado de Attenborough y Van Renterghem,
capítulo 4: las ecuaciones de gobierno (Ecs. 4.3-4.4), la actualización
leapfrog (Ecs. 4.11-4.12), la condición de estabilidad de Courant (Ecs.
4.13-4.14), el contorno rígido (Ec. 4.32) y el contorno de impedancia real
independiente de la frecuencia (Ecs. 4.33-4.35), expuestos a través de
`fdtd_simulation`/`FDTD2D` con las fuentes `GaussianPulse`, `CWSource` y
`SignalSource`, sondas de presión y una `obstacle_mask` rasterizada, con
`damping` como única pérdida volumétrica. La sección 3 deduce las cuatro
reglas de dimensionado que la regla de resolución no cubre (el ancho de banda
de la fuente, la duración de la ejecución, el grosor de la esponja en celdas,
y el transitorio y la ventana DFT de una medición en régimen estacionario), y
la sección 2 enumera las seis comprobaciones que dicen si una ejecución sirve.
Validado contra oráculos en forma cerrada: los modos normales de sala rígida
(Kuttruff cap. 3), los tiempos de llegada en campo libre y el decaimiento
cilíndrico en 2D, el eco de imagen de pared rígida, el coeficiente de
reflexión de impedancia y la relación de dispersión de la sección 6, con el
orden de convergencia medido coincidiendo con el diseño de segundo orden del
esquema. La cadena de campo cercano a lejano de la sección 4: captura de
fasores sobre contorno cerrado (`add_contour_probe`) y la integral de
Kirchhoff-Helmholtz 2D (`far_field_from_contour`, Williams cap. 8), validada
contra los campos exactos de las fuentes lineales monopolo y dipolo, la
propiedad de extinción, y los paneles QRD y metadifusor mallados frente a los
propios modelos de campo lejano Fraunhofer y TMM de la biblioteca.
La geometría interior es rígida por construcción: solo los cuatro lados del
dominio aceptan una impedancia o una esponja, así que una superficie interior
solo puede ablandarse respaldándola con una región `damping` con pérdidas, y
esa región es un fluido equivalente independiente de la frecuencia y no un
modelo poroso. El esquema es solo bidimensional: modela una fuente lineal en
sección transversal (propagación cilíndrica $1/\sqrt{r}$), no la propagación
esférica $1/r$ de una fuente puntual 3D, así que los niveles absolutos y las
tasas de decaimiento no son los de una sala 3D. El contorno abierto es la capa
absorbente de rampa cuadrática descrita como «el precursor sencillo» de una
verdadera capa perfectamente adaptada, no una PML en sí. Las ecuaciones de
gobierno suponen un medio en reposo, así que no se modela la advección por
viento o flujo, y el único contorno de impedancia es el real independiente de
la frecuencia de las Ecs. 4.33-4.35. Los medios elásticos no se cubren aquí:
el esquema compañero P-SV, sus superficies libres y su física fluido-sólido
tienen su propia guía, [Ondas elásticas y acoplamiento
fluido-sólido](/phonometry/es/simulation/elastic-waves/).
## Véase también
- [Ondas elásticas y acoplamiento fluido-sólido](/phonometry/es/simulation/elastic-waves/):
el esquema compañero P-SV sobre la misma malla escalonada, con ondas
de Rayleigh y de Scholte, conversión de modo y placas sumergidas.
- [Informe de conformidad](/phonometry/es/reference/conformance/): dos de los anclajes
de validación en forma cerrada de la sección 6 se ejecutan ahí.
- Referencia de la API: [`simulation.fdtd`](/phonometry/es/reference/api/simulation/fdtd/),
[`simulation.ntff`](/phonometry/es/reference/api/simulation/ntff/).
## Respuestas rápidas
### ¿Cómo elijo el paso de malla en FDTD?
Resuelve al menos 10 celdas por longitud de onda más corta,
$\Delta x \le c_\text{min} / (10 f_\text{max})$, con la velocidad de sonido
más baja del
dominio. Con exactamente 10 celdas la cota del error de dispersión sobre el
eje es alrededor del 1,6 %, reducida a en torno al 1,4 % con el
``cfl = 0.6`` por defecto (el número de Courant por eje es entonces
$S = \mathrm{cfl}/\sqrt{2} = 0{,}424$, así que el factor $1 - S^2$ vale 0,82), y
reducir $\Delta x$ a la mitad divide el error entre
cuatro (el esquema es de segundo orden). Con $\Delta x = 1$ cm el punto de 10
celdas en aire queda hacia los 3,4 kHz. Esa cota es un error de velocidad de
*fase*; un pulso viaja a la velocidad de grupo y va tres veces más lento
todavía, el 4,1 % con diez celdas, y lo paga por metro recorrido, que es lo
que [la sección 6 mide en tres mallas](#dispersion-numerica-y-precision).
### ¿Qué anchura deben tener el pulso de fuente y la esponja?
El espectro del `GaussianPulse` es $e^{-(\pi f\,\text{width})^2}$, 20 dB por
debajo en $f = 0{,}483/\text{width}$, así que con
`width` $\ge 0{,}5/f_\text{max}$ la energía inyectada se mantiene dentro de
la banda que la malla resuelve. La esponja es una regla en **celdas**, no en
longitudes de onda: como su tasa de absorción actúa a la vez sobre la
presión y sobre la velocidad, la capa conserva una impedancia real $\rho c$
y solo refleja su *gradiente*. Veinte celdas devuelven unos −61 dB a
incidencia normal, planos dentro de 3 dB de 125 Hz a 2 kHz; cuarenta celdas
devuelven unos −76 dB. Apretar `sponge_reflection` más allá de su valor por
defecto de 1e-4 hace la rampa más abrupta y empeora el residuo. La sección 3
desarrolla todo esto, y con ello la duración de la ejecución y la ventana
DFT.
### ¿Qué número de Courant mantiene estable una simulación FDTD?
El esquema explícito solo es estable mientras un frente de onda cruce como
mucho una celda por paso temporal. Con celdas cuadradas el número de Courant
es $\mathrm{CN} = c\,\Delta t\,\sqrt{2}/\Delta x \le 1$ (Attenborough y Van
Renterghem, Ec. 4.13). `fdtd_simulation` deriva el paso temporal del
parámetro `cfl` (ese mismo $\mathrm{CN}$, por defecto 0,6) y de la mayor
velocidad del sonido del mapa, y rechaza los valores fuera de $(0, 1)$. El
número de Courant *por eje* de la relación de dispersión de más arriba es un
número distinto, $S = c\,\Delta t/\Delta x = \mathrm{CN}/\sqrt{2}$, que vale
0,424 con el valor por defecto.
### ¿Puedo fiarme de los niveles absolutos de una ejecución FDTD 2D?
No. Una fuente puntual 2D es físicamente una fuente lineal infinita: su
amplitud se expande cilíndricamente como $1/\sqrt{r}$, 3,0 dB por
duplicación de distancia, en lugar del $1/r$ esférico (6,0 dB por
duplicación) de una fuente puntual 3D. Los patrones de interferencia y
difracción son fieles, pero los niveles absolutos y las tasas de caída no
son los de una sala 3D; valida cualquier afirmación cuantitativa 3D contra
una forma cerrada o un cálculo 3D.