はじめに
音響や電磁気を数値シミュレーションする方法の一つにFDTD法があります。これは偏微分方程式の微分を直接に差分に置き換えて計算するもので、理論もプログラムも簡単になる利点があります。
FDTD法は有限の領域を計算するので境界の取り扱いが重要です。
境界を0にしたままにすると固定端となり、波は全反射します。
開けた野外の音響など、無響室のように反射してほしくないときは境界に工夫がいります。
今回は吸収境界の中でもたぶん一番簡単なMurの一次吸収境界を導出し、実装して試してみます。
この分野を専門にしているわけではないので、間違っている可能性も十分あります。
その時はコメントで教えていただけるとありがたいです。
もとの画像はこんなふうに境界にノイズが乗っていないのですが、オンラインGIF圧縮ツールを使ったらこうなってしまいました。
参考
波動方程式の離散化はこちらを参考にしました。
Qiita:波動方程式の数値解法
Murの吸収境界はこちらを参考にしました。
光波の電磁界解析法の現状 時間領域差分法-光学分野への応用を期待して-
FDTD法
波動方程式の離散化
今回は二次元の波動方程式を考えます。
ここで$p$は圧力、$c$は音速です。
\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
中心差分で離散化します。
\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)
左辺に次の時刻の$p$がくるように式変形します。ここで$\Delta x$と$\Delta y$は等しいことにします。
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プログラム 固定端
プログラムの作成と実行はJupyter labで行いました。
import numpy as np
import matplotlib.pyplot as plt
from numba import jit
#単位はm,s,m/s
#dx = dyにする(式変形でそうしたから)
x_len = 8.0 #x方向の長さ m
y_len = 8.0 #y方向の長さ m
nx = 400 #x方向の格子数
ny = 400 #y方向の格子数
dx = x_len / nx #x方向の格子幅 m
dy = y_len / ny #y方向の格子幅 m
dt = 0.000001 #時間刻み s
c = 340.0 #音速 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("dxとdyを同じにする")
if c*dt/dx > 1.0:
print("CFL={} < 1".format(c*dt/dx))
def get_nearest(xp,yp):
"""
音源の置き場所計算用
"""
return np.argmin((x-xp)**2),np.argmin((y-yp)**2)
sx,sy = get_nearest(4.0,4.0) #音源の位置
@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):
#加振
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)) #現在のタイムステップ
p_pre = np.zeros((nx,ny)) #1つ前のタイムステップ
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
0.2msに1枚グラフを保存してGIMPでアニメーションにしました。
境界を0のままにすると、固定端反射になります。
Murの一次吸収境界 導出
一次元の波動方程式を考えます。
\frac{\partial^2 p}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0
$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
上式の第一項は後退波を表しているそうで、反射がないなら0になります。
\frac{\partial p}{\partial x} + \frac{1}{c} \frac{\partial p}{\partial t} = 0
これを差分化していきます。空間は$i-\frac{1}{2}$、時間は$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
$i-\frac{1}{2}$や$n+\frac{1}{2}$の位置に格子点はないので平均値を使います。
\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
あとは式変形して$p^{n+1}_i$を左辺に、それ以外を右辺に持っていきます。
まず両辺に$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
4種類の$p$があるので、それぞれまとめます。
\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
左辺の第一項以外を右辺に移して
\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}
両辺を$\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}
無事に参考文献と同じ形の式にできました。
FDTDプログラム Murの一次吸収境界
変更した部分だけ書きます。
def update(p,p_pre,t):
#加振
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の1次吸収境界
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
垂直に入射した波はよく消えますが、斜めに入射したものは少し残ってしまいます。
境界のノイズはGIF圧縮のせいです。
Murの二次吸収境界 導出
二次元の波動方程式に対して同様の操作を行えば、Murの二次吸収境界が導出できます。
\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
偏微分記号をたくさん書くのが大変なので、3つの項をそれぞれ$a$,$b$,$c$とおいて式変形します。
$a^2 + b^2 - c^2 = 0$の形を$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
$y$方向の後退波がない条件は
b + c \sqrt{1- \left( \frac{a}{c} \right) ^2} = 0
参考文献はここから一気に式変形が終わってしまうのですが、どうやら平方根の項をテイラー展開で近似しているようです。
$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
両辺に$c$をかけて
bc + c^2 - \frac{1}{2} a^2 = 0
記号をもとに戻して
\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
これがMurの二次吸収境界の条件式になります。
差分化も一次と同様にやろうと思ったら$p^{n+1}_{i-1,j}$など未来の位置の項が残ってしまい、うまくいきませんでした。しかも参考文献の差分化をそのままプログラムしたら発散するし、よくわかりません。
まとめ
FDTD法で重要な吸収境界のうちの一つを導出し実装しました。
回折などより精度が必要な現象を扱うときは、より反射が少ないPML吸収境界を使用するそうです。機会があればそちらも学んで投稿したいと思います。