[PYTHON] Simulation acoustique Méthode FDTD Dérivation et mise en œuvre des limites d'absorption de Mur

introduction

La méthode FDTD est l'une des méthodes de simulation numérique de l'acoustique et de l'électromagnétique. Ceci est calculé en remplaçant directement le différentiel de l'équation aux dérivées partielles par la différence et présente l'avantage de simplifier la théorie et la programmation.

Étant donné que la méthode FDTD calcule une aire finie, il est important de gérer la frontière. Si vous laissez la limite à 0, ce sera une extrémité fixe et l'onde sera totalement réfléchie. Si vous ne voulez pas qu'il soit réfléchi comme une pièce sans son, comme un son extérieur ouvert, vous devez définir une limite. Cette fois, nous allons dériver, implémenter et essayer la limite d'absorption primaire de Mur, qui est probablement la plus simple des limites d'absorption.

Je ne suis pas spécialisé dans ce domaine, il est donc fort possible que je me trompe. À ce moment-là, je vous serais reconnaissant si vous pouviez me le faire savoir dans les commentaires.

L'image d'origine n'a pas de bruit sur la bordure comme celui-ci, mais lorsque j'utilise l'outil de compression GIF en ligne, cela ressemble à ceci.

référence

Je me suis référé ici pour la dispersion de l'équation des vagues. Qiita: solution numérique de l'équation d'onde

Je me suis référé ici à la limite d'absorption de Mur.

[État actuel de la méthode d'analyse du champ électromagnétique pour les ondes lumineuses Application de la méthode de différence de région temporelle au champ optique-](https://annex.jsap.or.jp/photonics/kogaku/public/27-11-kaisetsu5. pdf)

Méthode FDTD comme méthode d'analyse du champ magnétique photoélectrique

Méthode FDTD

Discrimination de l'équation des vagues

Cette fois, nous considérerons une équation d'onde bidimensionnelle. Où $ p $ est la pression et $ c $ est la vitesse du son.

\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2}- \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0

Distinguer par différence de centre.

\frac{p^{n+1}_{i,j} - 2p^n_{i,j} + p^{n-1}_{i,j}}{\Delta t^2} = c^2 \left( \frac{p^{n}_{i+1,j} - 2p^n_{i,j} + p^{n}_{i-1,j}}{\Delta x^2} + \frac{p^{n}_{i,j+1} - 2p^n_{i,j} + p^{n}_{i,j-1}}{\Delta y^2} \right)

La formule est transformée de sorte que $ p $ à la prochaine fois vienne sur le côté gauche. Supposons maintenant que $ \ Delta x $ et $ \ Delta y $ soient égaux.

p^{n+1}_{i, j} = 2p^n_{i,j} - p^{n-1}_{i,j} + \frac{\Delta t^2 c^2}{\Delta x^2} \left( p^n_{i-1,j} + p^n_{i+1,j} + p^n_{i,j-1} + p^n_{i,j+1} - 4p^n_{i,j} \right)

Fin fixe du programme FDTD

Le programme a été créé et exécuté dans le laboratoire Jupyter.

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

#L'unité est m,s,m/s
#dx =Faire mourir(Parce que j'ai fait ça avec la transformation de formule)
x_len = 8.0 #Longueur dans la direction x m
y_len = 8.0 #Longueur dans la direction y m
nx = 400 #Nombre de grilles dans la direction x
ny = 400 #Nombre de grilles dans la direction y
dx = x_len / nx #Largeur du treillis dans la direction x m
dy = y_len / ny #Largeur du treillis dans la direction y m
dt = 0.000001 #Pas de temps s
c = 340.0 #Vitesse du son m/s
C1 = (dt**2)*(c**2)/(dx**2)
C2 = (c*dt)/dx
x = np.linspace(0.0,x_len,nx)
y = np.linspace(0.0,y_len,ny)
X,Y = np.meshgrid(x,y)

if not dx == dy:
    print("Faites dx et dy de la même manière")

if c*dt/dx > 1.0:
    print("CFL={} < 1".format(c*dt/dx))

def get_nearest(xp,yp):
    """
Pour calculer l'emplacement de la source sonore
    """
    return np.argmin((x-xp)**2),np.argmin((y-yp)**2)

sx,sy = get_nearest(4.0,4.0) #Position de la source sonore

@jit
def sound_source(p,p_pre,t,f):
    p[sx,sy] = np.sin(2.0*np.pi*f*t)
    p_pre[sx,sy] = np.sin(2.0*np.pi*f*(t-dt))

@jit
def update(p,p_pre,t):
    #Vibration
    f = 1.0e3
    if t < 1.0/f * 5.0:
        sound_source(p,p_pre,t,f)

    p_next = np.zeros_like(p)
    for i in range(1,nx-1):
        for j in range(1,ny-1):
            p_next[i,j]  = 2.0*p[i,j] - p_pre[i,j] + C1*(p[i-1,j] + p[i+1,j] + p[i,j-1] + p[i,j+1] - 4.0*p[i,j])
    return p_next,p

p = np.zeros((nx,ny)) #Pas de temps actuel
p_pre = np.zeros((nx,ny)) #Pas de temps précédent

t = 0
img_num = 1
for i in range(50000):
    if i%200==0:
        fig,axes = plt.subplots(figsize=(8,8))
        axes.imshow(np.flipud(p.T),vmin=-0.1,vmax=0.1,cmap="jet")
        fig.savefig("anime/{}.png ".format(img_num))
        plt.close()
        img_num += 1
    p,p_pre = update(p,p_pre,t)
    t += dt

J'ai enregistré un graphique toutes les 0,2 ms et je l'ai animé avec GIMP.

Laisser la limite à 0 entraîne des réflexions d'extrémité fixes.

Dérivation de la limite d'absorption primaire de Mur

Considérons une équation d'onde unidimensionnelle.

\frac{\partial^2 p}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0

Factoriser sous la forme $ a ^ 2-b ^ 2 = (a + b) (a-b) $.

\left( \frac{\partial p}{\partial x} + \frac{1}{c} \frac{\partial p}{\partial t} \right) \left( \frac{\partial p}{\partial x} - \frac{1}{c} \frac{\partial p}{\partial t} \right) = 0

Le premier terme de l'équation ci-dessus semble représenter une onde de recul, et il sera égal à 0 s'il n'y a pas de réflexion.

\frac{\partial p}{\partial x} + \frac{1}{c} \frac{\partial p}{\partial t} = 0

Nous allons différencier cela. L'espace est à $ i- \ frac {1} {2} $ et le temps est à $ n + \ frac {1} {2} $.

\frac{p^{n+\frac{1}{2}}_{i} - p^{n+\frac{1}{2}}_{i-1}}{\Delta x} + \frac{1}{c} \frac{p^{n+1}_{i-\frac{1}{2}} - p^n_{i-\frac{1}{2}}}{\Delta t} = 0

Puisqu'il n'y a pas de points de grille aux positions de $ i- \ frac {1} {2} $ et $ n + \ frac {1} {2} $, la valeur moyenne est utilisée.

\frac{ \frac{p^{n+1}_{i}+p^n_{i}}{2} - \frac{p^{n+1}_{i-1}+p^n_{i-1}}{2}}{\Delta x} + \frac{1}{c} \frac{\frac{p^{n+1}_{i} + p^{n+1}_{i-1}}{2} - \frac{ p^n_{i} + p^n_{i-1}}{2}}{\Delta t} = 0

Après cela, transformez la formule et amenez $ p ^ {n + 1} _i $ sur le côté gauche et les autres sur le côté droit.

Tout d'abord, multipliez les deux côtés par $ 2c \ Delta t \ Delta x $

c \Delta t \left( p^{n+1}_i + p^n_i - p^{n+1}_{i-1} - p^{n+1}_{i-1} \right) + \Delta x \left( p^{n+1}_i + p^{n+1}_{i-1} - p^n_i - p^n_{i-1} \right) = 0

Il existe quatre types de $ p $, je vais donc résumer chacun d'eux.

\left( c \Delta t + \Delta x \right) p^{n+1}_i
+ \left( c \Delta t - \Delta x \right) p^n_i
- \left( c \Delta t - \Delta x \right) p^{n+1}_{i-1}
- \left( c \Delta t + \Delta x \right) p^n_{i-1} = 0

Déplacer tout sauf le premier terme du côté gauche vers le côté droit

\left( c \Delta t + \Delta x \right) p^{n+1}_i = 
 \left( c \Delta t - \Delta x \right) \left( p^{n+1}_{i-1} - p^n_i \right)
 + \left( c \Delta t + \Delta x \right) p^n_{i-1}

Divisez les deux côtés par $ \ left (c \ Delta t + \ Delta x \ right) $

 p^{n+1}_i = 
 \frac{\left( c \Delta t - \Delta x \right)}{\left( c \Delta t + \Delta x \right)} \left( p^{n+1}_{i-1} - p^n_i \right)
 + p^n_{i-1}

La formule a la même forme que la référence.

Limite d'absorption primaire du programme FDTD Mur

N'écrivez que la pièce modifiée.

def update(p,p_pre,t):
    #Vibration
    f = 1.0e3
    if t < 1.0/f * 5.0:
        sound_source(p,p_pre,t,f)

    p_next = np.zeros_like(p)
    for i in range(1,nx-1):
        for j in range(1,ny-1):
            p_next[i,j]  = 2.0*p[i,j] - p_pre[i,j] + C1*(p[i-1,j] + p[i+1,j] + p[i,j-1] + p[i,j+1] - 4.0*p[i,j])

    #Limite d'absorption primaire de Mur
    for i in range(nx):
        p_next[i,0] = (c*dt-dx)/(c*dt+dx)*(p_next[i,1]-p[i,0]) + p[i,1]
        p_next[i,-1] = (c*dt-dx)/(c*dt+dx)*(p_next[i,-2]-p[i,-1]) + p[i,-2]

    for j in range(ny):
        p_next[0,j] = (c*dt-dx)/(c*dt+dx)*(p_next[1,j]-p[0,j]) + p[1,j]
        p_next[-1,j] = (c*dt-dx)/(c*dt+dx)*(p_next[-2,j]-p[-1,j]) + p[-2,j]

    return p_next,p

Les ondes incidentes verticalement disparaissent souvent, mais les ondes incidentes en diagonale restent un peu. Le bruit de bordure est dû à la compression GIF.

Dérivation de la limite d'absorption secondaire de Mur

La limite d'absorption quadratique de Mur peut être dérivée en effectuant la même opération sur l'équation d'onde bidimensionnelle.

\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2}- \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0

Puisqu'il est difficile d'écrire beaucoup de symboles différentiels partiels, la formule est transformée en définissant les trois termes comme respectivement $ a $, $ b $ et $ c $.

Changez la forme de $ a ^ 2 + b ^ 2-c ^ 2 = 0 $ en $ b ^ 2- \ left (c ^ 2-a ^ 2 \ right) = 0 $

b^2 - \left[  c^2 \left\{ 1 -  \left( \frac{a}{c} \right) ^2 \right\} \right] = 0

\left( b + c \sqrt{1- \left( \frac{a}{c} \right) ^2} \right) \left( b - c \sqrt{1- \left( \frac{a}{c} \right) ^2} \right) = 0

La condition qu'il n'y a pas de vague de recul dans la direction $ y $

b + c \sqrt{1- \left( \frac{a}{c} \right) ^2} = 0

Dans la référence, la transformation de formule se termine aussitôt à partir de là, mais il semble que le terme de la racine carrée soit approché par l'expansion de Taylor.

Quand $ x \ ll 1 $, $ \ sqrt {1-x ^ 2} \ approx 1- \ frac {x ^ 2} {2} $

b + c \left\{ 1 - \frac{1}{2} \left( \frac{a}{c} \right) ^2\right\} = 0

Multipliez $ c $ des deux côtés

bc + c^2 - \frac{1}{2} a^2 = 0

Remettez le symbole

\frac{\partial p}{\partial y \partial t} + \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} - \frac{c^2}{2} \frac{\partial^2 p}{\partial x^2} = 0

C'est l'expression conditionnelle de la limite d'absorption secondaire de Mur.

Quand j'ai essayé de faire la différenciation de la même manière que le primaire, cela n'a pas fonctionné parce que le terme de la position future telle que $ p ^ {n + 1} _ {i-1, j} $ est resté. De plus, si vous programmez la numérisation des références telle quelle, elle divergera et je ne suis pas sûr.

Résumé

Nous avons dérivé et mis en œuvre l'une des limites d'absorption importantes dans la méthode FDTD. Lorsqu'il s'agit de phénomènes nécessitant plus de précision, comme la diffraction, il semble utiliser la limite d'absorption PML avec moins de réflexion. Si j'ai une chance, j'aimerais apprendre et l'afficher.