[PYTHON] Akustische Simulation FDTD-Methode Ableitung und Implementierung der Mur-Absorptionsgrenze
Einführung
Die FDTD-Methode ist eine der Methoden zur numerischen Simulation von Akustik und Elektromagnetik. Dies wird berechnet, indem das Differential der partiellen Differentialgleichung direkt durch das Differential ersetzt wird, und hat den Vorteil, dass Theorie und Programmierung vereinfacht werden.
Da die FDTD-Methode eine endliche Fläche berechnet, ist es wichtig, die Grenze zu behandeln. Wenn Sie die Grenze bei 0 belassen, ist dies ein festes Ende und die Welle wird vollständig reflektiert. Wenn Sie nicht möchten, dass es wie ein nicht klingender Raum reflektiert wird, z. B. ein offener Außenklang, müssen Sie eine Grenze festlegen. Dieses Mal werden wir Mur's primäre Absorptionsgrenze ableiten, implementieren und ausprobieren, die wahrscheinlich die einfachste der Absorptionsgrenzen ist.
Ich bin nicht auf diesen Bereich spezialisiert, daher ist es gut möglich, dass ich falsch liege. Zu diesem Zeitpunkt würde ich mich freuen, wenn Sie mich in den Kommentaren wissen lassen könnten.
Das Originalbild weist am Rand kein Rauschen auf, aber wenn ich das Online-GIF-Komprimierungswerkzeug verwende, sieht es so aus.
Referenz
Ich habe hier auf die Streuung der Wellengleichung Bezug genommen. Qiita: Numerische Lösung der Wellengleichung
Ich bezog mich hier auf die Absorptionsgrenze von Mur.
[Aktueller Stand der Methode zur Analyse elektromagnetischer Felder für Lichtwellen Zeitbereichsdifferenzmethode - erwartete Anwendung auf das optische Feld](https://annex.jsap.or.jp/photonics/kogaku/public/27-11-kaisetsu5. pdf)
FDTD-Methode als photoelektrische Magnetfeldanalyse
FDTD-Methode
Diskriminierung der Wellengleichung
Dieses Mal werden wir eine zweidimensionale Wellengleichung betrachten. Wobei $ p $ Druck und $ c $ Schallgeschwindigkeit ist.
\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
Unterscheiden durch Mitteldifferenz.
\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)
Die Formel wird so transformiert, dass $ p $ beim nächsten Mal auf die linke Seite kommt. Nehmen wir nun an, dass $ \ Delta x $ und $ \ Delta y $ gleich sind.
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)
FDTD-Programm festes Ende
Das Programm wurde im Jupyter-Labor erstellt und ausgeführt.
import numpy as np
import matplotlib.pyplot as plt
from numba import jit
#Die Einheit ist m,s,m/s
#dx =Machen Sie dy(Weil ich das mit der Formeltransformation gemacht habe)
x_len = 8.0 #Länge in x-Richtung m
y_len = 8.0 #Länge in y-Richtung m
nx = 400 #Anzahl der Gitter in x-Richtung
ny = 400 #Anzahl der Gitter in y-Richtung
dx = x_len / nx #Gitterbreite in x-Richtung m
dy = y_len / ny #Gitterbreite in y-Richtung m
dt = 0.000001 #Zeitschritt s
c = 340.0 #Schallgeschwindigkeit 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("Machen Sie dx und dy das gleiche")
if c*dt/dx > 1.0:
print("CFL={} < 1".format(c*dt/dx))
def get_nearest(xp,yp):
"""
Zur Berechnung des Standorts der Schallquelle
"""
return np.argmin((x-xp)**2),np.argmin((y-yp)**2)
sx,sy = get_nearest(4.0,4.0) #Position der Schallquelle
@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)) #Aktueller Zeitschritt
p_pre = np.zeros((nx,ny)) #Vorheriger Zeitschritt
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
Ich habe alle 0,2 ms ein Diagramm gespeichert und es mit GIMP animiert.
Das Belassen der Grenze bei 0 führt zu Reflexionen am festen Ende.
Ableitung der primären Absorptionsgrenze von Mur
Betrachten Sie eine eindimensionale Wellengleichung.
\frac{\partial^2 p}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0
Faktor in die Form $ 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
Der erste Term der obigen Gleichung scheint eine zurückgehende Welle darzustellen, und er wird 0 sein, wenn keine Reflexion vorliegt.
\frac{\partial p}{\partial x} + \frac{1}{c} \frac{\partial p}{\partial t} = 0
Wir werden dies unterscheiden. Der Raum ist bei $ i- \ frac {1} {2} $ und die Zeit ist bei $ 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
Da an den Positionen von $ i- \ frac {1} {2} $ und $ n + \ frac {1} {2} $ keine Gitterpunkte vorhanden sind, wird der Durchschnittswert verwendet.
\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
Transformieren Sie danach die Formel und bringen Sie $ p ^ {n + 1} _i $ auf die linke Seite und die anderen auf die rechte Seite.
Multiplizieren Sie zunächst beide Seiten mit $ 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
Es gibt vier Arten von $ p $, daher werde ich jede zusammenfassen.
\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
Verschieben Sie alle bis auf den ersten Term auf der linken Seite nach rechts
\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}
Teilen Sie beide Seiten durch $ \ 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}
Die Formel hat die gleiche Form wie die Referenz.
FDTD-Programm Mur Primärabsorptionsgrenze
Schreiben Sie nur den geänderten Teil.
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])
#Mur's primäre Absorptionsgrenze
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
Vertikal einfallende Wellen verschwinden oft, aber diagonal einfallende Wellen bleiben ein wenig. Randrauschen ist auf GIF-Komprimierung zurückzuführen.
Ableitung der sekundären Absorptionsgrenze von Mur
Die quadratische Absorptionsgrenze von Mur kann abgeleitet werden, indem dieselbe Operation an der zweidimensionalen Wellengleichung durchgeführt wird.
\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
Da es schwierig ist, viele partielle Differentialsymbole zu schreiben, wird die Formel transformiert, indem die drei Terme als $ a $, $ b $ bzw. $ c $ festgelegt werden.
Ändern Sie die Form von $ a ^ 2 + b ^ 2-c ^ 2 = 0 $ in $ 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
Die Bedingung, dass es keine zurückgehende Welle in Richtung $ y $ gibt
b + c \sqrt{1- \left( \frac{a}{c} \right) ^2} = 0
In der Referenz endet die Formeltransformation sofort von hier aus, aber es scheint, dass der Term der Quadratwurzel durch Taylor-Expansion angenähert wird.
Wenn $ x \ ll 1 $ ist, ist es $ \ sqrt {1-x ^ 2} \ ca. 1- \ frac {x ^ 2} {2} $
b + c \left\{ 1 - \frac{1}{2} \left( \frac{a}{c} \right) ^2\right\} = 0
Multiplizieren Sie $ c $ auf beiden Seiten
bc + c^2 - \frac{1}{2} a^2 = 0
Setzen Sie das Symbol zurück
\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
Dies ist der bedingte Ausdruck für Mur's sekundäre Absorptionsgrenze.
Als ich versuchte, die Differenzierung auf die gleiche Weise wie die primäre durchzuführen, funktionierte dies nicht, da der Term der zukünftigen Position wie $ p ^ {n + 1} _ {i-1, j} $ erhalten blieb. Wenn Sie die Digitalisierung von Referenzen so programmieren, wie sie ist, wird sie außerdem abweichen, und ich bin mir nicht sicher.
Zusammenfassung
Wir haben eine der wichtigen Absorptionsgrenzen in der FDTD-Methode abgeleitet und implementiert. Bei Phänomenen, die mehr Präzision erfordern, wie z. B. Beugung, scheint die PML-Absorptionsgrenze mit weniger Reflexion verwendet zu werden. Wenn ich eine Chance habe, würde ich sie gerne lernen und posten.