Propagación de ondas acústicas por elementos finitos con R

La simulación por elementos finitos FDTD (Finite-Difference Time-Domain) intercambia la complejidad de las ecuaciones diferenciales que gobiernan determinados sistemas físicos (las ondas electromagnéticas, el calor, el sonido, los fluidos), por la fuerza bruta de los ordenadores. Se basa en discretizar en el espacio y en el tiempo las ecuaciones diferenciales asociadas a cada fenómemo, y en este artículo vamos a implementar las ecuaciones que modelan la propagación de ondas acústicas.

Al modelar la propagación del sonido en un medio (aire, agua,...), se está describiendo la propagación de pequeñas perturbaciones o vibraciones mecánicas locales, que sin embargo no conllevan desplazamiento neto de materia. Esas perturbaciones se pueden caracterizar principalmente mediante dos campos:

• La velocidad instantánea (magnitud vectorial) con la que las partículas del medio oscilan alrededor de su posición de equilibrio.
• La presión acústica (magnitud escalar) como variación de presión respecto a la presión estática del medio.

La presión es la más interesante a la hora de estudiar la propagación del sonido y será la única que abordemos en este ejercicio. A continuación se muestra de forma muy resumida la Ecuación de onda acústica referida a la presión, así como la forma de discretizarla de forma sencilla para el caso 2D que es el que vamos a estudiar (hacer clic para ver con más detalle):

Nuestro solver divide el espacio en una matriz con determinado número de celdas de ancho y de alto, y el tiempo también es discretizado en saltos parametrizables por el usuario. El valor de presión se actualiza en cada posición de la matriz en función de los valores de presión de esa misma celda y de celdas vecinas en el instante inmediatamente anterior y el previo a éste. Las ecuaciones se refieren a las posiciones por sus coordenadas (i,j) y a los instantes de tiempo por el valor de n (no es un exponente).

Las simulaciones por elementos finitos de fenómenos ondulatorios parecen mágicas. Con implementaciones aproximadas muy sencillas de ecuaciones complejas, podemos ver aparecer ante nuestros ojos todos los fenómenos aplicables a una onda: la propagación y su correspondiente atenuación, las interferencias (constructivas y destructivas), reflexiones, difracción, refracción, efecto doppler,...

Aprovechando que no hace mucho estudiamos en 'Simulando arrays con R' el comportamiento como array de agrupaciones lineales, vamos a simular la propagación de ondas y las distribuciones de lóbulos y nulos en arrays de altavoces isotrópicos (altavoces ideales que radian por igual en todas las direcciones), que podremos validar comparándolas con el diagrama polar del factor de array.

A continuación se muestran 5 simulaciones con un número creciente de elementos radiantes desde 1 (solo un altavoz) hasta 5. Los elementos siempre están separados λ/2 de sus vecinos, lo que ocasiona nulos en la dirección paralela al array cuando el número de altavoces es par, y pequeños lóbulos en esa dirección cuando este número es impar. Se ha eliminado el signo de las presiones ya que lo que nos interesa es ver el progreso de los frentes de onda así como las interferencias constructivas y destructivas entre contribuciones (hacer clic para ver en alta resolución):

Estas simulaciones son válidas para corroborar cualitativamente la física de la propagación del sonido, pero no son trasladables al mundo real por estar realizadas en su versión 2D. Eso significa que no siguen las atenuaciones que un sistema de audio va a sufrir en el mundo real, donde cualquier frente de ondas esférico se atenúa al propagarse según una proporcionalidad 1/r² (4πr² es la supeficie de una esfera), mientras que en nuestros modelos 2D la atenuación es por 1/r (longitud lineal 2πr).

Es muy interesante fijarse en las irregularidades que presentan los frentes de onda más externos (los primeros en generarse), y que se debena que en esos frentes no todos los elementos radiantes han participado por no haber tenido tiempo físico de llegar su contribución. Esas bandas anómalas donde la interferencia destructiva de los nulos no se cumple podemos considerarla un transitorio del cálculo.

Si uno mira el código verá que he tenido que añadir a la función generada y acelerada en C++ por ChatGPT, un offset de fase de π/2 en las fuentes de alimentación sinusoidales ya que de lo contrario se introducía un nivel de contínua DC de presión que llevaba a un desequilibrio entre las crestas positivas y negativas de los frentes de onda.

~~~

Repositorio con el código R: GitHub.