Simulación de ondas FDTD 2D

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 solucionador 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, promocionado a API pública con fuentes, sondas de presión, obstáculos rasterizados, condiciones de contorno por lado y un objeto de resultado congelado.

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 v (Attenborough y Van Renterghem, Ecs. 4.3-4.4):

$$\frac{\partial p}{\partial t}+\rho c^2\nabla \cdot \mathbf{v}=0,\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 (el número de Courant, por defecto 0,6) y de la mayor velocidad del sonido del mapa; los valores fuera de (0, 1) se rechazan porque el esquema carece de sentido más allá del límite (Ec. 4.14).

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, 0.01, 2.0e-3, shape=(200, 300),
    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 del §3)

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 senoidal 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.

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 difusores de cualquier forma son simples arrays booleanos. 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).
• 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 - ρc)/(Z + ρc); Z = ρc es anecoico.

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:

import numpy as np
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

3. Cuándo usarlo, y los límites del 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), donde un método en el dominio de la frecuencia necesita una resolución por frecuencia.

Mostrar el código de esta figura

import numpy as np
import matplotlib.pyplot as plt
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), la forma cerrada es miles de veces más barata: este solucionador es el contraste y el demostrador, no el sustituto. El oráculo funciona en ambos sentidos; una ejecución en caja rígida reproduce los modos analíticos de la sala:

import numpy as np
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]
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 exactamente 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 desplazamiento apenas visible hacia la izquierda de los picos más altos es la dispersión numérica del §4: las longitudes de onda cortas viajan ligeramente lentas en la malla, así que las resonancias modeladas leen por debajo en una fracción de porcentaje.

Mostrar el código de esta figura

import numpy as np
import matplotlib.pyplot as plt
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()

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. Trate las ejecuciones 2D como secciones y demostraciones, y valide cualquier afirmación cuantitativa 3D contra una forma cerrada o un solucionador 3D.

4. 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 dx)^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), con S = c dt / dx; 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 CN = 1. La regla práctica es resolver al menos 10 celdas por longitud de onda más corta, dx <= c_min / (10 f_max) con la velocidad de sonido más baja del dominio: con exactamente 10 celdas la cota de Courant pequeño (k dx)^2 / 24 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 dx = 1 cm el punto de 10 celdas en aire queda hacia los 3,4 kHz, y reducir dx 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). Los tests anclan el solucionador a oráculos analíticos: 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 anterior 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.

Qué cubre esta guía

Cubierto. 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. 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 del §4, con el orden de convergencia medido coincidiendo con el diseño de segundo orden del esquema.

No cubierto. El solver 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 asumen 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.

Véase también

• Informe de conformidad: dos de los anclajes de validación en forma cerrada del §4 se ejecutan ahí.
• Referencia de la API: simulation.fdtd.

Respuestas rápidas

¿Cómo elijo el paso de malla en FDTD?

Resuelva al menos 10 celdas por longitud de onda más corta, dx <= c_min / (10 f_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, y reducir dx a la mitad divide el error entre cuatro (el esquema es de segundo orden). Con dx = 1 cm el punto de 10 celdas en aire queda hacia los 3,4 kHz.

¿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 (el número de Courant, por defecto 0,6) y de la mayor velocidad del sonido del mapa, y rechaza los valores fuera de $(0,1)$.

¿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; valide cualquier afirmación cuantitativa 3D contra una forma cerrada o un solucionador 3D.

Referencias

• Attenborough, K. y Van Renterghem, T. (2021). Predicting outdoor sound (2.ª ed.). CRC Press. https://doi.org/10.1201/9780429470806Capítulo 4: el modelo FDTD presión-velocidad de referencia implementado aquí, desde las ecuaciones de gobierno (4.3-4.4) y la actualización leapfrog escalonada (4.11-4.12) hasta la condición de Courant (4.13-4.14), el análisis del error de fase (4.15) y las condiciones de contorno rígida y de impedancia finita (4.32-4.35). El módulo implementa este método numérico de libro, no una norma de medida, con este capítulo como formulación de referencia citable.
• Kuttruff, H. (2016). Room acoustics (6.ª ed.). CRC Press. https://doi.org/10.1201/9781315372150La sección 3.5 sitúa los métodos ondulatorios en el dominio del tiempo entre las aproximaciones numéricas a la ecuación de onda en recintos, y el capítulo 3 da los modos normales de sala rígida usados como oráculo analítico.