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 fuerza bruta computacional. Consiste en aproximar las ecuaciones asociadas a cada fenómemo discretizándolas en el espacio y en el tiempo. 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, sólidos,...), se está describiendo la propagación de pequeñas perturbaciones o vibraciones mecánicas locales, que sin embargo no conllevan desplazamiento neto de materia. Estas 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.

Fuente: University of Alberta

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 resumida la discretización de la Ecuación de onda acústica referida a la presión y 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 discretas (i,j) y a los instantes de tiempo por el valor también discreto de n (no es un exponente).

La función principal para actualizar el campo de presiones se la he pedido a ChatGPT, que además ha encapsulado el bucle masivo que recorre todas las posiciones de la matriz de presiones en una rutina compilada en C++ para lograr tiempos de ejecución razonables.

En todas las simulaciones FDTD se debe satisfacer que los parámetros de entrada de velocidad de propagación (c), resolución de la malla simulada (dx, dy) y paso temporal (dt), cumplan una condición de estabilidad denominada condición de Courant-Friedrichs-Lewy (CFL) o simplemente condición de Courant. Su interpretación física no puede ser más intuitiva, se trata de que en ningún momento nuestro algoritmo "corra" a mayor velocidad de la que supone avanzar como máximo una celda del espacio simulado por unidad de tiempo simulado, y que en el caso que nos ocupa se traduce como:

Las simulaciones por elementos finitos de fenómenos ondulatorios son mágicas. Con implementaciones aproximadas muy sencillas de ecuaciones originalmente complejas, podemos ver aparecer ante nuestros ojos todos los fenómenos aplicables a una onda: la propagación y su correspondiente atenuación, los fenómenos de interferencia constructiva y destructiva, las reflexiones y la absorción, la difracción o la refracción, el efecto Doppler,...

En 'Transferencia de calor por elementos finitos con R (I). Simulador' abordé en la era pre ChatGPT la transmisión del calor por diferencias finitas, y unos pocos años atrás, mi PFC 'Impulse antenna design using the FDTD' consistió en simular la propagación electromagnética (ecuaciones de Maxwell). Ya solo me queda en la lista de pendientes la reina de las simulaciones por elementos finitos: la mecánica de fluidos (ecuaciones de Navier-Stokes), de las que la propagación del sonido puede verse como un caso particular simplificado.

~~~

Aprovechando que no hace mucho estudiamos en 'Simulando arrays con R' el comportamiento como array de agrupaciones lineales de elementos radiantes, 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 al cancelarse dos a dos, y lóbulos secundarios en esa dirección cuando este número es impar al sobrevivir una contribución tras cancelarse las demás. 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 y el fenómeno de interferencia, pero no son trasladables al mundo real por estar realizadas en su versión 2D. Eso significa que no siguen las atenuaciones que van a darse en un sistema sonoro real, donde cualquier frente de ondas esférico pierde energía al propagarse según una proporcionalidad 1/r² (4πr² = superficie de la esfera), mientras que en nuestros modelos 2D realmente tenemos frentes de onda cilíndricos atenuándose por 1/r (2πr = longitud de la circunferencia).

Las irregularidades que presentan los frentes de onda más externos (los primeros en generarse), se deben a que no todos los elementos radiantes han participado en ellos al no haber tenido tiempo físico de llegar su contribución. Esas bandas anómalas donde la interferencia esperada no se cumple representan un transitorio del cálculo. Aquí podemos ver a modo de anillos de Saturno y para el ejemplo de 5 elementos radiantes, cómo pasamos por 1, 2, 3, 4 y hasta las 5 contribuciones finales, generándose nulos de presión cuando el número es par y presiones efectivas cuando es impar:

He tenido que añadir a la función generada por ChatGPT un offset de fase de π/2 en las fuentes de alimentación sinusoidales, convirtiendo el seno en un coseno. De lo contrario se introduce un nivel de continua (DC) de presión, que lleva a un desequilibrio entre las crestas positivas y negativas de los frentes de onda lo cual no es posible partiendo de un espacio en reposo acústico. Confieso que no termino de entender ni el origen del problema, ni por qué funciona mi solución puramente intuitiva (pasar de cuartos de ciclo seno ++-- a coseno +--+).

Terminamos con una animación del ejercicio con 4 elementos. No solo es vistosa e hipnótica, sino útil para poder comprobar cómo las interferencias que se esperan del array no se cumplen al 100% en campo cercano (zonas próximas a los elementos radiantes), pero conforme nos vamos al campo lejano, que es para el que se construye toda la teoría de arrays, ya sí aparecen los nulos y la señal donde esperamos que estén (hacer clic para ver la animación):

~~~

Repositorio con el código R: GitHub.