Comparative study of models of impedance boundary conditions in acoustic problems

Resumen

En este artículo se presentan distintas soluciones para la implementación numérica de condiciones de contorno de impedancia (reactancia local) en problemas acústicos. Para ello se analizan 2 tipos de ecuaciones: las ecuaciones de Euler y la ecuación de ondas, y se estudian diferentes soluciones para los contornos tanto en algoritmos de diferencias finitas en el dominio del tiempo (FDTD) como en algoritmos pseudo-espectrales en el dominio del tiempo (PSTD). El análisis de las distintas propuestas numéricas existentes en la literatura se realiza mediante exhaustivos experimentos numéricos que permiten estudiar el comportamiento absorbente de las distintas condiciones de contorno en función de la frecuencia y del ángulo de las ondas incidentes. Este novedoso estudio comparativo permite al ingenierio acústico escoger el modelo numérico que más se adapte a sus necesidades.

Abstract

In this paper, different implementations of numerical locally reacting boundary conditions are studied for acoustic problems. In this comparative study we analyze two types of equations, the Euler equations and the wave equation. We also analyze both finite-differences time-domain (FDTD) algorithms, and pseudo-spectral time domain (PSTD) numerical schemes. We compare different numerical implementations existing in the literature by means of exhaustive numerical experiments. These numerical experiments allow for the study of the absorbing properties of the different schemes as a function of the frequency and the angle of the incident sound waves. This novel comparative study will help the acoustic engineer in order to choose the proper numerical scheme for his/her simulations.

Palabras clave

Acústica de salas ; Diferencias finitas ; Condiciones de contorno

Keywords

Room Acoustics ; Finite Differences ; Boundary Conditions

1. Introducción

La distribución del sonido en una habitación es un fenómeno complejo que depende de la geometría del recinto y de las propiedades absorbentes de las paredes, del techo y del suelo (así como de los materiales que los recubren) [1] . Además, fenómenos tan típicamente acústicos como la difracción o las interferencias provocan que la distribución de las variables acústicas dependa fuertemente de la posición y del tiempo.

Así pues, los fenómenos acústicos en general son muy complejos y resulta extremadamente complicado encontrar soluciones analíticas para caracterizar el campo acústico. Es por ello que el uso de ordenadores para predecir el comportamiento del campo acústico se ha convertido en una herramienta fundamental.

En general, las simulaciones por ordenador en campos como la acústica arquitectónica, la aeroacústica o la acústica medioambiental se dividen en 2 grandes grupos [2] : los métodos geométricos y los métodos numéricos que solucionan la ecuación de ondas. El primer grupo comprende un conjunto de algoritmos basados en la hipótesis de que las longitudes de onda del sonido son significativamente más pequeñas que las dimensiones de los objetos presentes en el dominio de simulación. De entre estos métodos, los más populares son los métodos de trazado de rayos (ray-tracing en inglés) [3] , los conocidos como image-source methods[4] y los más recientes métodos llamados beam tracing methods[5] . Todos estos métodos presentan la desventaja de que no reproducen de forma natural fenómenos ondulatorios tan importantes como las interferencias o la difracción. Así pues, los resultados obtenidos con estos algoritmos son imprecisos en el rango de frecuencias bajas, proporcionando predicciones erróneas en este régimen.

Es por ello que los métodos numéricos de integración de las ecuaciones diferenciales que gobiernan la evolución del campo acústico son la alternativa preferida por investigadores e ingenieros. Desde el punto de vista matemático se trata de resolver un problema con condiciones iniciales y con condiciones de contorno. Para resolver este tipo de problemas hay diferentes alternativas: desde el uso de elementos finitos (FEM) [6]  and [7] o elementos de frontera (BEM) [8] , métodos que originalmente fueron desarrollados para tratar problemas en el dominio frecuencial, hasta los métodos en diferencias finitas en el dominio del tiempo (FDTD) [9] o los conocidos como digital waveguide mesh methods (DWM) [10] usados habitualmente para resolver problemas transitorios.

Todos estos métodos proporcionan diferentes ventajas y desventajas en términos computacionales dependiendo de su coste y complejidad. Sin embargo, la diferencia más acusada estriba en su rango de aplicabilidad: los métodos como FEM o BEM se centran en el análisis frecuencial y, por tanto, proporcionan resultados para situaciones estacionarias; por otro lado, los métodos como FDTD o DWM se usan principalmente para el cálculo de la respuesta impulsional, que es un tema central en acústica arquitectónica [1] y acústica ambiental [11] . Cabe decir, sin embargo, que en los últimos años han aparecido variantes de FEM y de BEM para tratar problemas acústicos en el dominio del tiempo [12] .

Recientemente han aparecido como alternativa unos nuevos algoritmos conocidos como métodos pseudo-espectrales en el dominio del tiempo (PSTD) [13] que usan transformadas de Fourier para calcular las derivadas espaciales. Los algoritmos PSTD hacen un uso exhaustivo de la fast fourier transform[14] por lo que, en algunas situaciones, son mucho más rápidos que los clásicos FDTD. Otra ventaja que tienen los métodos PSTD respecto a los FDTD es el hecho de que presentan una dispersión muy baja (en la práctica se consideran isótropos) [13] , [15]  and [16] . En los últimos años, los algoritmos PSTD se han usado en diferentes campos tales como la propagación de ondas de sonido [17]  and [18] , el modelado de transductores piezoeléctricos [19] o la simulación de dispositivos fotónicos [20] .

Para la aplicación de los algoritmos PSTD en situaciones de interés práctico es necesaria la formulación de modelos numéricos para las condiciones de contorno. A lo largo de la última década se ha prestado mucha atención a este problema en el campo de los algoritmos FDTD, y en particular a la obtención de modelos de impedancia local para caracterizar paredes reflectantes. Básicamente, este modelo asume que la presión acústica y la componente normal del la velocidad de partícula están linealmente relacionadas por una impedancia acústica Z[21] . Esta relación permite obtener expresiones analíticas del coeficiente de reflexión, R , obtenido en el régimen en el que una onda plana impacta con una pared de estas características. El desarrollo de condiciones de contorno de impedancia local para algoritmos PSTD es un tema muy reciente [18]  and [22] y que permitirá el uso de estos algoritmos en múltiples campos de investigación.

En el presente artículo se presenta por primera vez un estudio comparativo detallado de los modelos de impedancia local tanto para algoritmos FDTD como PSTD aplicados tanto a las ecuaciones de Euler como a la ecuación de ondas. El artículo está organizado como sigue: en la sección 2 se explica en detalle lo que significa el modelo de impedancia local para las condiciones de contorno; en la sección 3 se presentan los experimentos numéricos que servirán para cuantificar el comportamiento de los distintos modelos numéricos; en la sección 4 se analizarán diferentes algoritmos para las ecuaciones de Euler tanto en FDTD como en PSTD; en la sección 5 se realiza el mismo estudio para la ecuación de ondas; finalmente, la sección 6 está destinada a las conclusiones y a las líneas de investigación abiertas.

2. Modelo de impedancia local

En general existen muchas clases de materiales con propiedades acústicas diferentes. Existen algunos materiales porosos capaces de absorber mucha energía [23] ; también existen materiales en los que aparece el fenómeno de reflexión acústica caracterizados por su impedancia superficial[24] . Para este último tipo de materiales es razonable asumir que las partículas del material tienen un comportamiento local e independiente unas con otras.

Consideremos el caso de una superficie plana caracterizada por una impedancia de contorno específica. Si una onda plana impacta contra una superficie uniforme de extensión infinita, una parte de la energía acústica se reflejará en forma de otra onda plana, de amplitud y fase distinta con respecto de la onda incidente. Los cambios en la amplitud y en la fase que suceden durante la reflexión de una onda plana vienen definidos por el factor de reflexión complejo, R ,

$R=\vert \vert R\vert \vert e^{\left(\iota \phi \right)},$

( 1)

donde ||R || ≤ 1, ι es la unidad imaginaria y ϕ es la fase. El factor de reflexión viene fuertemente influenciado por las propiedades acústicas del material. Además, tanto su valor absoluto como su fase son dependientes de la frecuencia y dirección de la onda incidente. En otras palabras, la intensidad de una onda reflejada se ve disminuida un factor ||R ||² comparada con la onda incidente, perdiendo la fracción 1 − ||R ||² de energía. Esta cantidad es conocida como «coeficiente de absorción acústico», α ,

$\alpha =1-\vert \vert R\vert \vert ^{2}.$

( 2)

Para una pared de reflectividad cero, R  = 0, el coeficiente absorción toma su valor máximo 1. Estas paredes son conocidas como paredes totalmente absorbentes. Si R  = 1 (en fase), se conocen como paredes «duras» o «rígidas»; en el caso que R  = −1 (fase inversa), se habla de paredes «blandas». En ambos casos no hay absorción del sonido.

Como se ha mencionado antes, R depende fuertemente de las propiedades acústicas del material. Más concretamente, todas las propiedades acústicas vienen definidas por la impedancia acústica Z . La impedancia Z se define, en el caso de una onda plana, como el cociente entre la amplitud compleja y la componente normal de la velocidad de la partícula v . En general, esta impedancia presenta una respuesta particular dependiendo de la forma de la onda incidente, además de que puede ser fuertemente dependiente de la frecuencia. En este artículo se considerará el caso en el que Z es independiente de la frecuencia.

Para ser más concretos, definamos un dominio dos-dimensional V . La presión acústica, p (x , t ), y la velocidad de partícula, v (x , t ), dentro del dominio serán caracterizadas por la ecuación de ondas (o, equivalentemente, por las ecuaciones de Euler), excepto en esas posiciones, x , que estén localizadas en el contorno ∂V (véase fig. 1 ). La relación temporal entre la presión acústica y la velocidad en ∂V viene dada por

$p({\boldsymbol {\mbox{x}}},t)=Z{\boldsymbol {\mbox{v}}}({\boldsymbol {\mbox{x}}},t)\cdot {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}},$

( 3)

donde ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}}}$ es el vector normal de la pared [1] y Z es una constante positiva real (i.e. la impedancia). Por otro lado, en las superficies de impedancia local se asume que la ecuación de conservación lineal de masa se cumple en la dirección normal:

${\frac {\partial {\boldsymbol {\mbox{v}}}({\boldsymbol {\mbox{x}}},t)\cdot {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}}}{\partial t}}=$$-{\frac {1}{\rho }}({\boldsymbol {\nabla }}\cdot {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}})p({\boldsymbol {\mbox{x}}},t).$

( 4)

donde ∇ es el vector gradiente en coordenadas cartesianas.

Figura 1. Reflexión producida por una superficie de impedancia, ∂V , en el dominio 2D, V   . ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}}}$ es el vector normal asociado al contorno.

Si una onda plana que viaja hacia una superficie de impedancia local impacta contra ella con un ángulo de incidencia θ , se define el coeficiente de reflexión R como el cociente entre la onda reflejada y la incidente, p_ref y p_inc respectivamente. Este cociente está relacionado con la impedancia acústica Z a través de [21] :

$R(\theta )={\frac {Zcos\theta -\rho c}{Zcos\theta +\rho c}},$

( 5)

donde ρ es la densidad del aire, c representa la velocidad de propagación del sonido y θ es el ángulo representado en la fig 1 .

A pesar de que se puedan utilizar tanto las ecuaciones de Euler como la de ondas para caracterizar la propagación del sonido en el aire, la elección entre una u otra ecuación diferencial condicionará la forma de la ecuación numérica del modelo de impedancia. Por ejemplo, si se utilizan las ecuaciones de Euler, las condiciones de contorno se deducirán directamente de las ecuaciones (3) y (4) . Por el contrario, si se utilizan algoritmos basados en la ecuación de ondas para describir la propagación del sonido en el aire, se deben introducir unas variaciones a las ecuaciones (3) y (4) para que puedan ser aplicadas a estos algoritmos concretos. Por ejemplo, si se introduce la ecuación (3) en la ecuación de conservación de masa, ecuación (4) , se obtiene una expresion más apropiada para la ecuación de ondas,

${\frac {\partial p({\boldsymbol {\mbox{x}}},t)}{\partial t}}=$$-{\frac {Z}{\rho }}({\boldsymbol {\nabla }}\cdot {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}})p({\boldsymbol {\mbox{x}}},t).$

( 6)

Nótese que la ecuación (6) solo depende de la presión acústica, en lugar de la ecuación (4) que relaciona la presión con la componente normal de la velocidad.

3. Experimento numérico

En esta sección se define el experimento numérico que se utilizará para testear la precisión de varias condiciones de contorno numéricas en los métodos FDTD y PSTD desarrollados para simular el modelo de impedancia local. El diseño del experimento está inspirado en Kelloniemi et al. [25]  and [26] y está detallado en [18] . El experimento consiste en una malla dos-dimensional rectangular con una fuente localizada en x_s . Varios receptores, ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ y ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }^{'}}}$ , donde ξ  = 1, 2, 3…, se sitúan a lo largo de líneas paralelas τ y τ ′ tal y como se muestra en la figura 2 .

Figura 2. Una representación esquemática del experimento numérico: la fuente se localiza en xs , los receptores se situan en τ y τ ′. ∂V representa una superficie plana que está a la misma distancia de xτ que de ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }^{'}}}$ .

Utilizando el experimento ilustrado en la figura 2 , 2 simulaciones distintas se llevan a cabo: una simulación incluyendo una línea de nodos de contorno, ∂V , localizado justo en el medio; y una segunda simulación en el espacio libre sin ninguna pared en medio. En ambas simulaciones, la fuente de sonido puntual, x_s , emite un impulso acústico que se aproxima por:

$p_{fuente}^{n}={\frac {sin(2\pi f(n-n_{0})\Delta t)}{2\pi f(n-n_{0})\Delta t}}\quad n\leq n_{t}.$

( 7)

Nótese que esta función tiene un espectro plano desde 0 hasta f . n representa el paso de tiempo y n_t es el tiempo en el que el impulso acústico cesa.

En la primera simulación, la presión acústica se mide en todas las posiciones, ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ , ξ  = 1, 2, 3…. Estas señales contienen no solo el sonido directo, sino también el sonido reflejado por la superficie, ∂V   . Por otro lado, en la segunda simulación (en el espacio libre) las señales se miden en ambas localizacions, en ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ y en ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }^{'}}}$ . De esta forma, se puede eliminar el sonido directo en los receptores ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ de la primera simulación utilizando la información obtenida en las mismas posiciones pero en la segunda simulación.

Finalmente, para poder medir el coeficiente de reflexión en los receptores ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ , la respuesta al impulso en el espacio de frecuencias obtenida en la primera simulación (una vez eliminado el sonido directo) se compara con los resultados obtenidos en la segunda simulación en las posiciones ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }^{'}}}$ .

En resumen: para cada receptor ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ es posible obtener el coeficiente de reflexión comparando el espectro de las señales medidas en las posiciones ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ con las posiciones especulares ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }^{'}}}$ . Debido al hecho de que distintos receptores corresponden a distintos ángulos, con este experimento somos capaces de medir en una única simulación el error absoluto

$\epsilon =20{log}_{10}\parallel R_{teo}-R_{med}\parallel ,$

( 8)

para cada frecuencia y ángulo, simplemente comparando el coeficiente de reflexión medido R_med con su valor teórico R_teo deducido de la ecuación (5) . Por último, la parte derecha de la función Hann se ha usado para promediar la mitad de la señal y así eliminar errores de truncamiento en el cálculo del espectro.

A lo largo del presente artículo consideraremos un dominio dos-dimensional V de 2000 × 2000 nodos con una pared de impedancia local localizada en los nodos (i , 1000). En este caso, la fuente de sonido se localizará en el nodo (520, 900). Además, la señal de entrada es un impulso acústico generado con una fuente blanda. Todos los casos han sido testeados con f  = 2.500 Hz, n₀  = 40, n_t  = 80, Δt  = 1/16.000 s en la ecuación (7) . En todos los casos se ha considerado el mínimo número de estabilidad de Courant permitido para cada esquema numérico, ya sea FDTD o PSTD. Los ángulos de incidencia considerados en el presente experimento están comprendidos en θ  ∈ [0, 80^∘ ], aproximadamente. Cabe mencionar que no hay una distribución homogénea de ángulos en este experimento pero que, por motivos ilustrativos, los ángulos desconocidos se han obtenido mediante interpolaciones lineales. Finalmente, para evitar errores de contorno que pudiesen afectar a la medida, la simulación se ha llevado a cabo en hasta 1.024 iteraciones temporales.

4. Condiciones de contorno para las ecuaciones de Euler

En esta sección consideraremos las ecuaciones de Euler en 2 dimensiones, usadas en multitud de problemas acústicos ya que describen fielmente la propagación espacio-temporal del campo acústico [21] en entornos cerrados:

${\frac {\partial v_{x}}{\partial t}}+{\frac {1}{\rho }}{\frac {\partial p}{\partial x}}=0$ ${\frac {\partial v_{y}}{\partial t}}+{\frac {1}{\rho }}{\frac {\partial p}{\partial y}}=0$ ${\frac {\partial p}{\partial t}}+\rho c^{2}\left[{\frac {\partial v_{x}}{\partial x}}+{\frac {\partial v_{y}}{\partial y}}\right]=0,$

( 9)

donde ρ es la densidad del aire y c es la velocidad de propagación del sonido. A partir de ahora, consideraremos que tanto ρ como c   son constantes. La velocidad de las partículas del aire viene dada por el vector ${\textstyle {\boldsymbol {\mbox{v}}}=}$$(v_{x},v_{y})$ y p es la presión acústica relativa. Las primeras 2 ecuaciones afirman que un gradiente de presión produce una aceleración del fluido (en este caso, el aire), mientras que la tercera ecuación nos dice que la divergencia de la velocidad produce una compresión en el fluido. Finalmente, decir que estas ecuaciones son válidas para pequeñas velocidades y para valores pequeños de la presión relativa.

4.1. Modelos de impedancia local para algoritmos FDTD

En 1995 [9] se presentó una expresión numérica de condiciones generales de contorno para el algoritmo leap-frog formulado en una malla escalonada. Los algoritmos basados en este tipo de mallas tienen la característica de que tanto la presión acústica como la velocidad de partícula se calculan alternativamente en posiciones y tiempos discretos distintos. Es decir, la presión p (x , y , t   ) y las 2 componentes de la velocidad ${\textstyle v_{x}(x,y,t)}$ y ${\textstyle v_{y}(x,y,t)}$ pasan a ser las cantidades discretas ${\textstyle p\vert _{i,j}^{n}}$ , ${\textstyle v_{x}\vert _{i+1/2,j}^{n+1/2}}$ y ${\textstyle v_{y}\vert _{i+j+1/2}^{n+1/2}}$ . Consecuentemente, por como está definido el algoritmo, en la superficie ∂V (ver fig. 2 ) solo podrá haber nodos que computen ${\textstyle v_{x}\vert _{i+1/2,j}^{n+1/2}}$ .

Las ecuaciones numéricas que describen el comportamiento de impedancia local se derivan de las ecuaciones (3) y (4) presentadas en la sección 2 . Para obtener el esquema numérico se asume un operador asimétrico de diferencias finitas para las derivadas espaciales. El resultado final es una ecuación capaz de simular tanto impedancias que dependen de la frecuencia como impedancias constantes. Este documento se focaliza en el estudio de modelos de impedancia local Z   independientes de la frecuencia. El esquema numérico para los nodos de velocidad ${\textstyle v_{x}\vert _{i+1/2,j}^{n+1/2}}$ que pertenecen a ∂V tiene la siguiente forma

$v_{x}\vert _{i+1/2,j}^{n+1/2}=\gamma v_{x}\vert _{i+1/2,j}^{n-1/2}-$$\beta p\vert _{i,j}^{n},$

( 10)

con

$\gamma ={\frac {1-Z/Z_{FDTD}}{1+Z/Z_{FDTD}}}$ $\beta ={\frac {1}{1+Z/Z_{FDTD}}}$ $Z_{FDTD}={\frac {\rho \Delta }{\Delta t}}.$

( 11)

Notar que γ y β son constantes adimensionales que dependen muy fuertemente de la impedancia de los materiales.

Debido al hecho de que Z  ∈ [0, + ∞], es más conveniente mostrar los resultados utilizando el coeficiente de reflexión en la dirección normal, R_n . Es decir, la ecuación (5) cuando θ  = 0 ya que R_n está definido entre [− 1, 1]. Para ver el comportamiento de las condicones de contorno, se han estudiado los casos R_n  = −1, − 0,9, …, 1 (ΔR_n  = 0,1). Finalmente, se ha medido el error absoluto comparando los resultados experimentales con el coeficiente de reflexión teórico, R_teo , expresado en dB.

Las simulaciones para medir el coeficiente de reflexión se han llevado a cabo con Δt  = 1/16.000 y el número de stabilidad de Courant óptimo para el algoritmo leap-frog   , ${\textstyle S=}$$1/{\sqrt {2}}$ . Estos resultados se muestran en la figura 3 . Cada gráfica representa el error absoluto fijado por la ecuación (8) en función del ángulo de incidencia, θ , y la frecuencia. Este error se representa en una escala de grises donde el negro representa valores de pocos dB negativos y el color blanco equivale a errores menores de −40 dB.

Figura 3. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (10) , combinado con el algoritmo leap-frog . De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

Cada gráfica viene etiquetada por el coeficiente de reflexión normal que, como se ha dicho antes, fija un valor de la impedancia Z a través de la ecuación (5) con θ  = 0. Por tanto, para cada valor de R_n (i.e. de Z   ) se ha obtenido el error absoluto en distintos ángulos, correspondientes a los receptores ${\textstyle {\boldsymbol {\mbox{x}}}_{{\tau }_{\xi }}}$ , para todas las frecuencias menores de 2.500 Hz.

Como se puede ver, la precisión de este modelo de impedancia local es muy alta observando errores menores de 30 dB en casi todo el rango de R_n y para todas las frecuencias. Solo para R_n  → 0,3, este error aumenta de forma homogénea a −20 dB, que es un error suficientemente pequeño para ser ignorado. En conclusión, este modelo es apropiado para simular el modelo de impedancia local en todo el rango de Z . A pesar de todo, se debe mencionar que este algoritmo está definido en una malla escalonada, limitando considerablemente la definición de paredes, sobre todo en entornos complejos.

4.2. Modelo de impedancia local para el método pseudo-espectrales

La primera implemetación de condiciones de contorno parcialmente absorbentes para métodos pseudo-espectrales con las ecuaciones de Euler ha sido propuesta recientemente [22] . Para poder definir estas condiciones de contorno, se presenta un esquema PSTD basado en una malla centrada, lo que implica que las cantidades acústicas se evalúan en el mismo instante y en el misma posición. La idea de definir un algoritmo formulado en una malla centrada permite, de manera muy simple, crear un modelo de impedancia local simplemente teniendo en cuenta la dirección normal de la pared a caracterizar mediante un parámetro adimensional ξ . Así pues, las condiciones de evolución para aquellos nodos de la red que están en la frontera toman la siguiente forma:

• Para ξ ≤ 1 :

$v_{x}\vert _{i,j}^{n+1}=v_{x}\vert _{i,j}^{n}-{\frac {\Delta t}{\rho }}F_{x}^{-1}\left[\iota {\frac {2\pi n_{x}}{N_{x}\Delta }}F_{x}[p\vert _{:,j}^{n}]\right],$ $v_{y}\vert _{i,j}^{n+1}=v_{y}\vert _{i,j}^{n}-{\frac {\Delta t}{\rho }}F_{y}^{-1}\left[\iota {\frac {2\pi n_{y}}{N_{y}\Delta }}F_{y}[p\vert _{i,:}^{n}]\right],$ $p\vert _{i,j}^{n+1}=\xi \left(p\vert _{i,j}^{n}-\rho c^{2}\Delta tF_{x}^{-1}\left[\iota {\frac {2\pi n_{x}}{N_{x}\Delta }}F_{x}[v_{x}\vert _{:,j}^{n+1}]\right]\right),$

• Para ξ > 1 :

$v_{x}\vert _{i,j}^{n+1}={\frac {1}{\xi }}\left(v_{x}\vert _{i,j}^{n}-{\frac {\Delta t}{\rho }}F_{x}^{-1}\left[\iota {\frac {2\pi n_{x}}{N_{x}\Delta }}F_{x}[p\vert _{:,j}^{n}]\right]\right),$ $v_{y}\vert _{i,j}^{n+1}=v_{y}\vert _{i,j}^{n}-{\frac {\Delta t}{\rho }}F_{y}^{-1}\left[\iota {\frac {2\pi n_{y}}{N_{y}\Delta }}F_{y}[p\vert _{i,:}^{n}]\right],$ $p\vert _{i,j}^{n+1}=p\vert _{i,j}^{n}-\rho c^{2}\Delta tF_{x}^{-1}\left[\iota {\frac {2\pi n_{x}}{N_{x}\Delta }}F_{x}[v_{x}\vert _{:,j}^{n+1}]\right]$ $.$

( 12)

donde (i , j ) son las coordenadas espaciales; Δ  = Δx  = Δy   representa la discretización espacial de la malla dos-dimensional. ${\textstyle F_{\mu }}$ y ${\textstyle F_{\mu }^{-1}}$ denotan la transformada discreta de Fourier sobre el eje μ y su inversa respectivamente; n_x y n_y son los índices de las transformadas de Fourier; y el símbolo: denota todas las coordenadas según el eje μ   . Finalmente, ${\textstyle \iota =}$${\sqrt {-1}}$ y N_μ es el número total de puntos de la malla sobre el eje μ . La constante de estabilidad de Courant para este esquema numérico viene dada por [13] :

$S=c{\frac {\Delta t}{\Delta }}\leq {\frac {2}{\pi {\sqrt {D}}}},$

( 13)

donde D representa la dimensión.

De este esquema se observa que las 3 cantidades acústicas, p y v , están evaluadas en el mismo tiempo, n , y lugar (i , j ). Observamos que se mantiene el esquema PSTD en la frontera donde las derivadas espaciales se calculan mediante tranformadas discretas de Fourier (para más detalles de esta formulación, el lector puede remitirse a la referencia [22] ). Además, nótese que las ecuaciones que actualizan la presión acústica solo calculan el gradiente en la dirección x , ya que para este experimento concreto, la pared de impedancia local es paralela al eje y . Por último, se debe resaltar que el parámetro, ξ , está relacionado con el coeficiente de reflexión, R , a través de una función que depende fuertemente del número de estabilidad de Courant S . Esta relación ha sido obtenida a numéricamente [22] y su demostración analítica sigue, a día de hoy, abierta,

• Para ξ ≤ 1 :

$\xi ={\frac {1+R}{1+R+Scos(\theta )(1-R)-S(1+R)}},$

( 14)

• Para ξ > 1 :

$\xi ={\frac {RS+S+cos(\theta )(RS-S+1-R)}{cos(\theta )(1-R)}}.$

( 15)

Para verificar numéricamente la relación entre el parámetro ξ y el coeficiente de reflexión, ecuaciones (14) y (15) , se realizaron los experimentos explicados en la sección 3 . Para estos experimentos se ha utilizado un esquema PSTD formulado en una malla centrada combinado con las ecuaciones que caracterizan este modelo semiempírico [22] . Comparando el coeficiente de reflexión analítico obtenido a través de las ecuaciones (14) y (15) con los resultados del experimento, se obtendrá el error absoluto que genera el modelo para todos los ángulos y hasta 2.500 Hz. En la figura 4 se grafican los resultados para nueve valores de ξ : 0, 0,25, 0,5, 0,75, 1, 2,5, 5, 7,5 y 10 (de a ) a i ) respectivamente, cubriendo así todo el rango de valores de R_n .

Figura 4. La función error para distintos valores de ξ en función del ángulo de incidencia y la frecuencia. Las gráficas corresponden a: a )ξ  = 0; b )ξ  = 0,25; c )ξ  = 0,5; d )ξ  = 0,75; e )ξ  = 1; f )ξ  = 2,5; g )ξ  = 5; h )ξ  = 7,5; i )ξ  = 10.

Las simulaciones dan resultados muy buenos en todo el rango de ξ para ángulos de incidencia pequeños (θ  ≤ 40^∘ ) donde la desviación respecto el comportamiento analítico es menor de −30 dB. Las desviaciones más relevantes respecto al comportamiento analítico aparecen en las zonas en que ξ  → 1. A pesar de todo, el comportamiento en casi todo el rango de frecuencias y ángulos es aceptable. Por último, notar que para ξ → ∞, aparece un error a altas frecuencias que parece incrementar con el valor de ξ . Como veremos a continuación este error aparece única y exclusivamente por formular el esquema en PSTD.

4.3. Reseñas adicionales

En la sección anterior se acaba de presentar un modelo semi-empírico para el método de PSTD, definiendo un parámetro adimensional ξ y computando únicamente el gradiente en la dirección normal en la ecuación que actualiza la presión. Además, se ha presentado una expresión analítica que relaciona este parámetro con el coeficiente de reflexión, el ángulo y el número de estabilidad de Courant (ver ecuaciones (14) y (15) ). Estas relaciones se puede escribir en términos de la impedancia de la pared:

• Para ξ ≤ 1 :

$\xi ={\frac {Z/(\rho c)}{S+Z/(\rho c)-ZS/(\rho c)}},$

( 16)

• Para ξ > 1 :

$\xi =ZS/(\rho c)-S+1,$

( 17)

El comportamiento de reactancia local de las ecuaciones (16) y (17) se ha testeado en 2 dimensiones a través del experimento definido en la sección 3 . Los bajos errores mostrados por el método sugieren que este esquema puede ser utilizado en muchas situaciones prácticas de problemas acústicos con el método de PSTD. Cabe mencionar, sin embargo, que a pesar de los buenos resultados, este método debe ser complementado con condiciones de contorno perfectamente absorbentes, tales como las bien conocidas perfectly matched layers (PML) [27]  and [28] , con el fin de evitar el fenómeno de Gibbs que aparece cuando se deriva espectralmente una función periódica no continua [29] .

Por otro lado, otro factor muy importante a considerar de este esquema es su remarcable flexibilidad debido al hecho de que puede ser extendido de forma muy simple y directa a cualquier algoritmo euleriano FDTD formulado en una malla centrada. Por ejemplo, para obtener el esquema numérico de condiciones de contorno para el algoritmo leap-frog se debe sustituir los términos que computan la derivada espacial de la ecuación (12) por operadores de diferencias finitas backward/forward. Por tanto, el esquema numérico de condiciones de contorno para el método de FDTD estándar viene dado por,

• Para ξ ≤ 1 :

$v_{x}\vert _{i,j}^{n+1}=v_{x}\vert _{i,j}^{n}-a_{1}(p\vert _{i+1,j}^{n}-p\vert _{i,j}^{n}),$ $v_{y}\vert _{i,j}^{n+1}=v_{y}\vert _{i,j}^{n}-a_{1}(p\vert _{i,j+1}^{n}-p\vert _{i,j}^{n}),$ $p\vert _{i,j}^{n+1}=\xi [p\vert _{i,j}^{n}-a_{2}(v_{x}\vert _{i,j}^{n+1}-v_{x}\vert _{i-1,j}^{n+1})],$

• Para ξ > 1 :

$v_{x}\vert _{i,j}^{n+1}={\frac {1}{\xi }}[v_{x}\vert _{i,j}^{n}-a_{1}(p\vert _{i+1,j}^{n}-p\vert _{i,j}^{n})],$ $v_{y}\vert _{i,j}^{n+1}=v_{y}\vert _{i,j}^{n}-a_{1}(p\vert _{i,j+1}^{n}-p\vert _{i,j}^{n}),$ $p\vert _{i,j}^{n+1}=p\vert _{i,j}^{n}-a_{2}(v_{x}\vert _{i,j}^{n+1}-v_{x}\vert _{i-1,j}^{n+1}),$

( 18)

donde ${\textstyle a_{1}={\frac {\Delta t}{\rho \Delta }}}$ y ${\textstyle a_{2}=}$${\frac {\rho c^{2}\Delta t}{\Delta }}$ . Se puede observar que las derivadas espaciales de estas condiciones de contorno se calculan con operadores de diferencias finitas en lugar de las derivadas espectrales que se utilizaban en el esquema PSTD, ecuación (12) . Soprendentemente, el parámetro ξ que aparece en las condiciones de contorno (18) mantiene exactamente la misma relación con la impedancia acústica Z que la que mantenía el modelo formulado en PSTD (ver ecuaciones (16) y (17) ). Esto implica que estas condiciones de contorno son completamente independientes del método numérico utilizado para aproximar la derivada espacial de la EDP. Este hecho enfatiza aún más la tremenda utilidad que tendría encontrar una derivación analítica para demostrar la generalidad de las ecuaciones (16) y (17) .

Igual que antes, se ha realizado el mismo experimento para testear las propiedades absorbentes del esquema numérico. En este caso, la ecuación 18 se ha combinado con un algoritmo FDTD centrado. Como en el experimento presentado en la sección 4.1 , se ha fijado Δt  = 1/16.000 s y ${\textstyle S=}$$1/{\sqrt {2}}$ , que es el número de estabilidad máximo permitido para este método. Los resultados del experimento se muestran en la figura 5 . Se puede observar que el error absoluto es menor de −20 dB en casi todo el rango de ξ . Solo para ξ  = 1, se observan errores mayores de −20 dB a grandes ángulos, θ  > 40. Como en el modelo PSTD, el método presenta dificultades para simular zonas muy absorbentes (ver fig. 5e ). Al parecer, este error que aparece en las zonas que ξ  → 1 es totalmente independiente del método empleado e intrínseco del esquema semi-empírico. Por otro lado, los resultados (fig. 5 ) son claramente mejores en el rango de ξ → ∞ si se comparan con los obtenidos con el método PSTD (ver fig. 4 ).

Figura 5. La función error para distintos valores de ξ en función del ángulo de incidencia y la frecuencia. Las gráficas corresponden a: a )ξ  = 0; b )ξ  = 0,25; c )ξ  = 0,5; d )ξ  = 0,75; e )ξ  = 1; f )ξ  = 2,5; g )ξ  = 5; h )ξ  = 7,5; i )ξ  = 10.

No obstante, hay que mencionar que los resultados obtenidos con los métodos FDTD y PSTD no han sido realizados bajo las mismas condiciones. Básicamente, la diferencia radica en el número de estabilidad de Courant, que para FDTD es ${\textstyle S=}$$1/{\sqrt {2}}$ y para PSTD es ${\textstyle S=}$$2/(\pi {\sqrt {2}})$ . Esto implica que la discretización espacial, Δ , en FDTD se ve reducida un factor 2/π   , lo que lleva a mallados más refinados que en las simulaciones PSTD. En otras palabras, si las simulaciones FDTD se hubiesen realizado con un ${\textstyle S=}$$2/(\pi {\sqrt {2}})$ , se hubiera visto reflejado en un incremento de los errores absolutos en las mismas zonas que en los resultados PSTD.

5. Condiciones de contorno para la ecuación de ondas

En esta sección nos centraremos en el comportamiento de distintos modelos numéricos de impedancia local para la ecuación de ondas para la presión acústica:

${\frac {{\partial }^{2}p({\boldsymbol {\mbox{x}}},t)}{\partial t^{2}}}-$$c^{2}\vartriangle p({\boldsymbol {\mbox{x}}},t)=0,$

( 19)

También partiremos de la ecuación de reactancia local, ecuación (6) . La ecuación de ondas es la más utilizada en problemas de acústica de salas [1] y en simulaciones de acústica virtual [7] .

5.1. Modelo de impedancia local para algoritmos FDTD

La primera aproximación de condiciones de contorno en acústica de salas para el método de diferencias finitas fue presentada por Huopaniemi et al. [30] , siendo la primera definición 1-D para el método de guía de onda digital. Básicamente, este método define su ecuación numérica teniendo en cuenta el concepto de coeficiente de reflexión. El esquema numérico tiene la siguiente forma explicita

$p\vert _{i,j}^{n+1}=(1+R_{n})p\vert _{i-1,j}^{n}-R_{n}p\vert _{i,j}^{n-1}.$

( 20)

Se ha elegido estudiar estas condiciones de contorno porque el análisis de su comportamiento de reactancia local nunca se había llevado a cabo en una simulación 2-D. En este caso, los nodos que no pertenezcan a ∂V son actualizados con el algoritmo discreto FDTD basado en la ecuación de ondas fijando Δt  = 1/16.000 s y el número máximo de estabilidad S permitido por el algoritmo.

Los resultados se muestran en la figura 6 . Se puede observar que solo para R_n  ≥ 0,8 y R_n  ≤ −0,6 las condiciones de contorno numéricas, ecuación (20) , dan errores absolutos menores que −20 dB. Para el resto de R_n , el error absoluto es considerablemente mayor de −20 dB. Cabe mencionar que los materiales absorbentes son poco comunes en simulaciones de este tipo. En otras palabras, el coeficiente de absorción α (α  = 1 − ||R_n ||² ) varía solo entre 0 y 0,5. En conclusión, bajo estas circunstancias se puede afirmar que este modelo numérico tiene resultados relativamente aceptables para aplicaciones como acústica, acústica de salas, aeroacústica …

Figura 6. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (20) , combinado con un esquema FDTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

No obstante, existen otras condiciones de contorno presentadas por Kowalczyk y van Walstijn [31] que mejoran incluso las condiciones de contorno presentadas en la figura 3 . Este modelo está basado en 2 ecuaciones: la ecuación de ondas y la ecuación (6) . Ambas ecuaciones se aproximan mediante operadores de diferencias centradas. La ecuación resultante se puede escribir en términos de un nodo definido fuera de V , conocido como ghost point , que para una pared de las características de nuestro experimento viene dada por

$p\vert _{i,j}^{n+1}={\frac {1}{\left(1+Z\lambda \right)}}2(1-2{\lambda }^{2})p\vert _{i-1,j}^{n}+{\lambda }^{2}(p\vert _{i,j+1}^{n-1}$ $+p\vert _{i,j-1}^{n-1})+2{\lambda }^{2}p\vert _{i-1,j}^{n}+(Z\lambda -1)p\vert _{i,j}^{n-1}.$

( 21)

Como se puede ver, esta ecuación depende fuertemente de un parámetro λ que para simulaciones 2D su número óptimo coincide con S   (i.e. ${\textstyle \lambda =}$$1/{\sqrt {2}}$ ). Luego se verán las limitaciones producidas por λ en simulaciones híbridas con algoritmos PSTD para los nodos de propagación (véase sección 5.3 ).

Los experimentos numéricos realizados con la ecuación (21) se han llevado a cabo bajo las mismas condiciones que los hechos para la ecuación (20) . Los resultados se muestran en la figura 7 , donde las gráficas de color muestran precisiones casi perfectas en todo el rango acústico de impedancia incluso para altas frecuencias. Por tanto, claramente estas condiciones son mucho mas adecuadas que las de la ecuación (20) , pudiendo simular incluso R_n  → 0 a muy alta precisión.

Figura 7. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (21) , combinado con un esquema FDTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

5.2. Modelo de impedancia local para algoritmos PSTD

En esta sección se presenta la única formulación realizada por el momento de un modelo de impedancia local para algoritmos PSTD [18] . La formulación del esquema numérico se basa en mezclar la técnica de PSTD para el algoritmo de propagación con una ecuación numérica construida con operadores de diferencias finitas para modelar el contorno. Como se muestra en la referencia [18] este algoritmo híbrido permite tener un esquema computacional capaz de modelar la propagación del sonido en espacios cerrados con un nivel de precisión bastante aceptable.

El punto de partida del esquema de condiciones de contorno de impedancia local es la ecuación (6) presentada en la sección 2 . Por motivos de estabilidad derivamos esta ecuación respecto del tiempo, con lo que nos queda:

${\frac {{\partial }^{2}p({\boldsymbol {\mbox{x}}},t)}{\partial t^{2}}}=$$-{\frac {Z}{\rho }}({\boldsymbol {\nabla }}\cdot {\overset {\mbox{ˆ}}{\boldsymbol {\mbox{n}}}}){\frac {\partial p({\boldsymbol {\mbox{x}}},t)}{\partial t}}.$

( 22)

Una vez alterada la ecuación (6) , el esquema numérico se obtiene utilizando operadores de diferencias finitas para aproximar los operadores de derivadas parciales de la ecuación (22) .

Para las derivadas temporales se utilizan los operadores de primera y segunda derivada centrados de segundo orden, y para la derivada espacial se utiliza un operador forward /backward dependiendo de la dirección normal de la pared a simular. En el caso concreto del experimento que se está tratando, el esquema numérico de impedancia local se expresa de la siguiente manera

$p\vert _{i,j}^{n+1}={\frac {2\rho \Delta }{\rho \Delta +0,5Z\Delta t}}p\vert _{i,j}^{n}-{\frac {\rho \Delta -0,5Z\Delta t}{\rho \Delta +0,5Z\Delta t}}p\vert _{i,j}^{n-1}$ $+{\frac {0,5\Delta tZ}{\rho \Delta +0,5Z\Delta t}}(p\vert _{i-1,j}^{n+1}-p\vert _{i-1,j}^{n-1}).$

( 23)

A diferencia de otros modelos basados en la ecuación (6) , este modelo es estable en todo el rango Z . Esta conclusión se saca tras hacer un análisis de Von Neumann que también se puede encontrar en [18] .

Con el objetivo de estudiar la idoneidad de las condiciones de contorno, ecuación (23) , combinado con una simulación multidimensional de PSTD, se han llevado a cabo distintos experimentos numéricos de acuerdo con lo definido en la sección 3 . El tiempo de discretización se ha fijado a Δt  = 1/16.000 s y el número de estabilidad de Courant ha sido el máximo permitido para el método de PSTD cuyo valor es ${\textstyle S=}$$2/(\pi {\sqrt {2}})$ .

La concordancia de los resultados con las predicciones teóricas es bastante buena en la mayoría de casos tratados. Del rango R_n  = −1 a R_n  = −0,3 se pueden considerar resultados muy buenos, ya que el error absoluto es ϵ  ≤ −25 dB para cualquier ángulo y frecuencia. Para el rango R_n  = −0,2 a R_n  = 0,2 el error absoluto crece a altas frecuencias y ángulos de incidencia pequeños. Por otro lado, de R_n  = 0,2 a R_n  = 0,8 el error absoluto vuelve a ser más que aceptable ya que alcanza valores ϵ  ≤ −20 dB. Finalmente, para R_n  = 0,9 y R_n  = 1 el error absoluto es homogéneo con errores del orden de −15 dB, los cuales se consideran demasiado grandes para ser un error aceptable.

Estos resultados no son para nada inesperados: por un lado, las regiones con valores del error absoluto mayores que −10 dB coinciden con aquellas regiones cuyo coeficiente de reflexión teórico, R_teo , resulta menor que −35 dB en escala logarítmica. Es sabido que estas regiones casi absorbentes son muy difíciles de caracterizar mediante cualquier método numérico multidimensional (se puede encontrar en la literatura técnica por ejemplo en [15] ); por otro lado, el aumento del error absoluto en el rango de R_n  > 0,8 sugiere que el uso de algoritmos híbridos (PSTD para la propagación y FDTD para el modelo de impedancia local) introduce un error numérico inherente que llega a ser inaceptable solo para R_n  → 1.

Para poder entender mejor el comportamiento de las condiciones de contorno (23) , se ha realizado el mismo experimento pero utilizando el algoritmo de propagación FDTD. Las simulaciones se han llevado a cabo con el mismo Δt   , pero con el número de estabilidad distinto y consecuentemente con distinta discretización espacial. Para estos experimentos, se ha fijado el número de estabilidad de Courant óptimo para el algoritmo FDTD, que es ${\textstyle S=}$$1/{\sqrt {2}}$ .

Los resultados se presentan en la figura 9 . Se observa que la precisión mejora considerablemente respecto de los resultados obtenidos con el algoritmo híbrido. Más concretamente, los resultados son claramente mejores en el rango de R_n  > 0,3. Por tanto, se puede afirmar que la combinación de FDTD para el contorno y PSTD para el algoritmo de propagación introduce un error inherente que, como hemos visto, para R_n  → 1 es suficientemente crítico para que sea considerado inaceptable. Por otro lado, el error absoluto otra vez es mayor cuando el coeficiente de reflexión tiende a 0, al contrario de propuestas como la de [31] que dan resultados casi perfectos en todo el rango de R_n , incluso para valores muy absorbentes. Por tanto, se puede concluir que, independientemente del algoritmo de propagación, las condiciones de contorno, ecuación (23) , dan sus mayores errores cuando R_n es cercano a 0. El error que aparece en el rango R_n  > 0,3 es debido a la formulación de un esquema numérico híbrido.

Figura 9. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (23) combinada con la ecuación de ondas discreta en FDTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

5.3. Discusión

Como se acaba de ver, las condiciones de contorno, ecuación (23) , son apropiadas para ser utilizadas en simulaciones PSTD. A pesar de todo, también se ha visto que debido al hecho de que el esquema tiene una formulación híbrida, aparece un error que a partir de un cierto rango de R_n incrementa y que es solamente crítico en R_n  = 1. Además, la ecuación (23) presenta dificultades intrínsecas para simular zonas más absorbentes. Por tanto, una posible duda que aparece sería saber si la ecuación (23) es la mejor opción en simulaciones PSTD, teniendo en cuenta que los resultados obtenidos con la formulación entera en FDTD también presentaban un incremento del error en las zonas más absorbentes (ver fig. 9 ) que no aparecían en los presentados en la sección 5.1 (ver fig. 7 ).

Llegados a este punto y siguiendo con la misma estrategia de formulación híbrida, sería interesante saber si las formulaciones presentadas en la sección 5.1 combinadas con un algoritmo PSTD son más apropiadas que las presentadas por Spa et al. [18] . Primeramente, se analizan las condiciones de contorno ecuación (20) combinadas con un algoritmo PSTD que, como en el experimento de la sección anterior, Δt  = 1/16.000 s y ${\textstyle S=}$$2/\pi {\sqrt {2}}$ . A pesar de los malos resultados obtenidos con este modelo en un esquema FDTD (ver fig. 6 ) y sumado al error que se espera por formular un modelo híbrido, creemos oportuno estudiar su comportamiento para sacar conclusiones más variadas sobre las formulaciones híbridas en aplicaciones acústicas.

Los resultados se presentan en la figura 10 . Como se esperaba, el error absoluto en casi todo el rango de R_n es mayor de −15 dB, el cual es demasiado alto como para ser ignorado. Cabe mencionar que el error se incrementa en todo R_n respecto a los resultados obtenidos con el esquema puramente FDTD. Por tanto, esto corrobora lo dicho antes de que las formulaciones híbridas introducen un error numérico que en este caso concreto es crítico en todo el rango de R_n .

Figura 10. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (20) combinada con la ecuación de ondas discreta en PSTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

Con el objetivo de mejorar los resultados obtenidos con la ecuación (20) , se han testeado las condiciones de contorno, ecuación (21) , propuestas por Kovalczyk y van Walstijn [31] en una simulación PSTD. Por un lado, se esperan mejores resultados que los obtenidos con la ecuación (20) ya que los resultados obtenidos con este esquema combinados con un algoritmo FDTD han tenido los resultados más precisos de este artículo. Por otro lado, un problema que se deriva de la formulación ecuación (21) es que la elección del parámetro λ no es unívoca debido a la fuerte relación entre el parámetro λ y el número de estabilidad de Courant del algoritmo de propagación PSTD.

Por tanto, se han llevado a cabo 2 experimentos combinando la ecuación (21) con un algoritmo PSTD para caracterizar la propagación. En una simulación, se ha fijado ${\textstyle \lambda =}$$1/{\sqrt {2}}$ y en la otra, ${\textstyle \lambda =}$$2/(\pi {\sqrt {2}})$ . La elección de los coeficientes está directamente relacionada con el número S de cada formulación. Las figuras 11 y 12 ilustran los resultados obtenidos en los 2 experimentos. Sorprendentemente, en ambos casos, el error absoluto es bastante mayor en casi todo el rango de R_n , si se comparan con los resultados obtenidos con la ecuación (23) (ver fig. 8 ). La simulación que utiliza ${\textstyle \lambda =}$$1/{\sqrt {2}}$ obtiene mejores resultados que la que utiliza el número óptimo de PSTD. A pesar de todo, para el mejor caso, en términos de resultados, solo se tienen resultados aceptables en el rango de R_n  < −0,3 donde el error decrece de forma homogénea de −20 a −40 dB. Estos errores son bajos y comparables a los errores obtenidos con la ecuación (23) . No obstante, de estos resultados se concluye que, por el momento, las condiciones híbridas, ecuación 23 , son la mejor opción para simular problemas de acústica para PSTD.

Figura 11. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (21) , fijando ${\textstyle \lambda =}$$1/{\sqrt {2}}$ combinada con la ecuación de ondas discreta en PSTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

Figura 12. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (21) , fijando ${\textstyle \lambda =}$$2/(\pi {\sqrt {2}})$ combinada con la ecuación de ondas discreta en PSTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

Figura 8. El error absoluto expresado en decibelios para distintos valores del coeficiente de reflexión normal Rn obtenidos con el esquema numérico, ecuación (23) combinada con la ecuación de ondas discreta en PSTD. De arriba abajo y de izquierda a derecha: a) Rn  = 1, b) Rn  = 0,9, c)Rn  = 0,8, d) Rn  = 0,7, e) Rn  = 0,6, f) Rn  = 0,5, g) Rn  = 0,4, h) Rn  = 0,3, i) Rn  = 0,2, j) Rn  = 0,1, k) Rn  = 0, m) Rn  = −0,1, n) Rn  = −0,2, o) Rn  = −0,3, p) Rn  = −0,4, q) Rn  = −0,5, r) Rn  = −0,6, s) Rn  = −0,7, t) Rn  = −0,8, u) Rn  = −0,9 y v) Rn  = −1.

6. Conclusiones

En el presente artículo se ha hecho un estudio comparativo de diferentes condiciones de contorno de impedancia local para problemas acústicos. Para hacer un análisis exhaustivo de las propiedades de absorción de los distintos modelos y compararlos entre sí, se han elaborado tests numéricos donde se calcula numéricamente el coeficiente de absorción para distintos ángulos y distintas frecuencias. El estudio se ha hecho para gran parte de los algoritmos existentes en la literatura tanto para la ecuación de ondas como para las ecuaciones de Euler. Las principales conclusiones de este estudio se pueden resumir:

• Para algoritmos FDTD en ecuaciones de Euler, las condiciones de contorno presentadas en referencia [9] se muestran eficaces y, por el momento, no hay en la literatura condiciones que las mejoren significativamente.
• Para algoritmos PSTD en ecuaciones de Euler solo hay un modelo de impedancia local en la literatura, referencia [22] , que proporciona resultados más que aceptables.
• Para algoritmos FDTD en la ecuación de ondas, las condiciones de contorno presentadas en la referencia [31] son insuperables en simulaciones en 2 dimensiones. Cabe decir que en 3 dimensiones las cosas son bastante distintas, ya que hay una ambigüedad en la implementación de este método. Además, hay que tener en cuenta que las simulaciones en 2 dimensiones son físicamente distintas de las simulaciones en 3 dimensiones [32] .
• Para algoritmos PSTD en la ecuación de ondas, el modelo híbrido FDTD-PSTD presentado en [18] ha demostrado ser el más eficaz.

Así pues, el presente trabajo representa una completa guía de modelos de condiciones de contorno para el científico o el ingeniero que quiera hacer simulaciones de problemas que se describen o bien con las ecuaciones de Euler o bien con la ecuación de ondas. Finalmente, enfatizar el hecho de que hacen falta nuevos modelos para algoritmos PSTD que mejores sus prestaciones.

Agradecimientos

El trabajo de A.G. está financiado en parte por el Ministerio de Industria a través del proyecto TSI-020100-2008-462 del Plan Avanza I+D. El trabajo de C.S. está parcialmente financiado por la agencia Chilena CONICYT dentro del proyecto FONDECYT num. 3110046. El trabajo de L.P. está parcialmente financiado por la agencia Chilena CONICYT dentro del proyecto FONDECYT num. 11100253.

References

1. [1] H. Kutruff; Room Acoustics; (4th Edition)Spon Press (2000)
2. [2] Savioja L. (1999). Modeling Techniques for Virtual Acoustics. PhD Thesis, Helsinki University of Technology, Helsinki, Finlandia.
3. [3] A. Krokstad, S. Strom, S. Sorsdal; Calculating the acoustical room response by the use of a ray tracing technique; J. Sound Vibration, 8 (1968), pp. 118–125
4. [4] J.B. Allen, D.A. Berkley; Image method for efficiently simulating small-room acoustics; J. Acoust. Soc. Am., 65 (1979), pp. 943–950
5. [5] T. Funkhouser, N. Tsingos, I. Carlbom, G. Elko, M. Sondhi, J. West, et al.; A beam tracing method for interactive architectural acoustics; J. Acoust. Soc. Am., 115 (2004), pp. 739–756
6. [6] J.R. Wright; An exact model of acoustic radiation in enclosed spaces; J. Audio Eng. Soc., 43 (1995), pp. 813–820
7. [7] L. Savioja, A. Järvinen, K. Melkas, K. Saarinen; Determination of the low frequency behaviour of an IEC listening room.; Proc. of the Nordic Acoustical Meeting (1996), pp. 55–58
8. [8] R.D. Ciskowski, C.A. Brebbia; Boundary element methods in acoustics; Computational Mechanics Publications Southampton Boston Co-published with Elsevier Applied Science (1991)
9. [9] D. Botteldooren; Finite-difference time-domain simulation of low-frequency room acoustic problems; J. Acoust. Soc. Am., 98 (1995), pp. 3302–3308
10. [10] D.T. Murphy, A. Kelloniemi, J. Mullen, S. Shelley; Acoustic modeling using the digital waveguide mesh; IEEE Signal Process. Mag., 24 (2007), pp. 55–66
11. [11] T.V. Renterghem, E.M. Salomons, D. Botteldooren; Efficient FDTD-PE model for sound propagation in situations with complex obstacles and wind profiles; Acta Acust. United Acust., 91 (2005), pp. 671–679
12. [12] J.A. Hargreaves, T.J. Cox; A transient boundary element method model of Schroeder diffuser scattering using well mouth impedance; J. Acoust. Soc. Am., 124 (5) (2008), pp. 2942–2951
13. [13] Q.H. Liu; The PSTD algorithms: a time-domain method requiring only two cells per wavelenght; Microw. Opt. Technol. Lett., 15 (1997), pp. 158–165
14. [14] J.W. Cooley, J.W. Tukey; An algorithm for the machine calculation of complex Fourier series; Math. Comput., 19 (1965), pp. 297–301
15. [15] A. Taflove, S.C. Hagness; Computational Electrodynamics: The Finite-Difference Time-Domain Method; Artech House Publishers (2000)
16. [16] C. Spa, T. Mateos, A. Garriga; Methodology for studying the numerical speed of sound in finite differences schemes; Acta Acust. United Acust., 95 (4) (2009), pp. 690–695
17. [17] Q.H. Liu; The pseudospectral time-domain (PSTD) algorithm for acoustic waves in absorptive media; IEEE trans. Ultrason. Ferroelectr. Freq. Control, 45 (4) (1998), pp. 1044–1055
18. [18] C. Spa, A. Garriga, J. Escolano; Impedance boundary conditions for pseudo-spectral time-domain methods in room acoustics; Appl. Acoust., 71 (5) (2010), pp. 402–410
19. [19] E. Filoux, S. Callé, D. Certon, M. Lethiecq, F. Levassort; Modeling of piezoelectric transducers with combined pseudospectral and finite-difference methods?; J. Acoust. Soc. Am., 123 (6) (2008), pp. 4165–4173
20. [20] W.H.P. Pernice; Pseudo-spectral time-domain simulation of the transmission and the group delay of photonic devices; Opt. Quant. Electron, 40 (2008), pp. 1–12
21. [21] P.M. Morse, K.U. Ingard; Theoretical Acoustics; Princeton University Press (1986)
22. [22] C. Spa, J. Escolano, A. Garriga; Semi-empirical boundary conditions for the linearized acoustic euler equations using pseudo-spectral time-domain methods; Appl. Acoust., 72 (4) (2011), pp. 226–230
23. [23] L.L. Beranek; Acoustic impedance of commercial materials and the perfomance of rectangular rooms with one treated surface; J. Acoust. Soc. Am., 12 (1940), pp. 14–27
24. [24] T. Embleton, J. Piercy, N. Olson; Outdoor sound propagation over ground of finite impedance; J. Acoust. Soc. Am., 59 (2) (1976), pp. 267–277
25. [25] A. Kelloniemi, L. Savioja, V. Välimäki; Spatial filter-based absorbing boundary for the 2-D digital waveguide mesh; IEEE Signal Process. Lett., 12 (2) (2005), pp. 126–129
26. [26] A. Kelloniemi; Improved adjustable boundary condition for the 2-D digital waveguide mesh; Proceedings of the 8th Int. Conference on Digital Audio Effects (DAFx 05), Madrid (2005), p. 2005
27. [27] J.-P. Berenger; Three-dimensional perfectly matched layer for the absorption of electromagnetic waves; J. Comput. Phys., 127 (1996), pp. 363–379
28. [28] X. Yuan, D. Borup, J.W. Wiskin, M. Berggren, R. Eidens, S. Johnson; Formulation and validation of Berengers PML absorbing boundary for the FDTD simulation of acoustic scattering; IEEE Trans. on Ultrasonics, Ferroelectrics and Frequency Control, 44 (1997), pp. 816–822
29. [29] B. Fornberg; A Practical Guide to Pseudospectral Methods; Cambridge University Press, Cambridge, UK (1996)
30. [30] J. Huopaniemi, L. Savioja, M. Karjalainen; Modeling of reflections and air absorption in acoustical spaces: a digital filter design approach; IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA’97), New Platz, NY, USA (1997)
31. [31] K. Kowalczyk, M. van Walstijn; Multichannel sound reproduction based on wave field synthesis; International Symposium on Room Acoustics. Satellite Symposium of the 19th International Congress on Acoustics, Seville, 2007 (2007)
32. [32] C. Spa, J. Escolano, A. Garriga, T. Mateos; Compensation of the afterglow phenomenon in 2-D discrete-time simulations; IEEE Signal Process. Lett., 17 (8) (2010), pp. 758–761