---
title: "Simulación de ondas FDTD 2D"
description: "Un solucionador acústico 2D determinista por diferencias finitas en el dominio del tiempo sobre una malla escalonada presión-velocidad: fuentes gaussianas, tonales y de señal arbitraria, sondas de presión, obstáculos rasterizados, contornos rígidos, de impedancia y absorbentes, y un resultado congelado con historias de sonda, instantáneas del campo y un método plot, siguiendo el capítulo 4 de Attenborough y Van Renterghem (2021)."
references:
  - type: book
    authors: ["Attenborough, K.", "Van Renterghem, T."]
    year: 2021
    title: "Predicting outdoor sound"
    edition: "2.ª ed."
    publisher: "CRC Press"
    doi: "10.1201/9780429470806"
    note: "Capí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."
  - type: book
    authors: ["Kuttruff, H."]
    year: 2016
    title: "Room acoustics"
    edition: "6.ª ed."
    publisher: "CRC Press"
    doi: "10.1201/9781315372150"
    note: "La 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."
---

import ThemeImage from '../../../../components/ThemeImage.astro';

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.

<ThemeImage src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/diagram_fdtd_es.svg" alt="Flujo desde la definición del dominio (mapas de velocidad del sonido y densidad con el paso de malla dx) y la geometría (máscara de obstáculos y condiciones de contorno por lado), pasando por las fuentes inyectadas en celdas, la actualización leapfrog en malla escalonada de velocidad y presión y la condición de estabilidad de Courant, hasta el FDTDResult congelado con historias de sonda, instantáneas del campo y un método plot" width="86%" loading="eager" />

## 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,
\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` (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).

```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, 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)
```

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:

```python
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.

<ThemeImage src="https://raw.githubusercontent.com/jmrplens/phonometry/main/.github/images/fdtd_simulation_es.webp" alt="Dos paneles. Izquierda: una instantánea del campo de presión en un dominio de 3 por 2 metros con una barrera rígida vertical delgada, mostrando el frente directo, la reflexión que vuelve hacia la fuente y la onda difractada en el borde de la barrera, con la fuente marcada con una estrella y dos sondas con puntos. Derecha: la historia de presión en ambas sondas; la sonda con visión directa muestra el pulso directo y la reflexión de la barrera, la sonda en sombra una llegada difractada más débil y retardada" width="96%" />

<details>
<summary>Mostrar el código de esta figura</summary>

```python
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)
res.plot(ax=ax_p)
plt.tight_layout()
plt.show()
```

</details>

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:

```python
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 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.

## Véase también

- [Informe de conformidad](/phonometry/es/reference/conformance/): dos de los anclajes
  de validación en forma cerrada del §4 se ejecutan ahí.
- Referencia de la API: [`simulation.fdtd`](/phonometry/es/reference/api/simulation/fdtd/).
