Mur 흡수 경계 조건
概要
波動の伝播を数値的に計算する方法には有限要素法(FDM)、$k$-space methodなどがある。これらの方法は波動が無限に伝播することを仮定しているが、実際のシミュレーションでは有限の格子内で波動が伝播する。そのため、境界では波動が反射する問題が生じる。これを解決するためには、実際にシミュレーションしたい領域よりもはるかに広く格子を設定するか、境界条件をうまく設定して波動の反射が起こらないようにする。そのような境界条件を吸収境界条件absorbing boundary condition (ABC)という。
公式
次のような波動方程式を考えよう。
$$ \dfrac{\partial ^{2}}{\partial t^{2}} u(x,t) = c^{2} \dfrac{\partial ^{2}}{\partial x^{2}} u(x,t) \quad \text{for } (x, t) \in \mathbb{R} \times [0, \infty) $$
これを有限な格子空間$[0, 1] \times [0, T]$でシミュレーションするとしよう。$x$軸空間を$M$等分すると$\Delta x = 1/M$、時間を$N$等分すると$\Delta t = T/N$である。波動方程式の解$u$の離散化を$U_{i}^{n} = u(i \ast \Delta x, n \ast \Delta t)$という。このとき$i = 0, 1, \cdots, M$、$n = 0, 1, \cdots, N$である。吸収境界条件は次のようになる。
$$ U_{0}^{n+1} = U_{1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{1}^{n+1} - U_{0}^{n} \right) $$
$$ U_{M}^{n+1} = U_{M-1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{M-1}^{n+1} - U_{M}^{n} \right) $$
説明
1次元と2次元で吸収境界条件を適用した結果は次の通りである。
導出1
波動方程式を次のように整理しよう。
$$ \begin{align*} && \dfrac{\partial ^{2}}{\partial t^{2}} u(x,t) &= c^{2} \dfrac{\partial ^{2}}{\partial x^{2}} u(x,t) \\ \implies && \dfrac{\partial ^{2}}{\partial t^{2}} u(x,t) - c^{2} \dfrac{\partial ^{2}}{\partial x^{2}} u(x,t) &= 0 \\ \implies && \left(\dfrac{\partial }{\partial t} - c \dfrac{\partial }{\partial x} \right) \left(\dfrac{\partial }{\partial t} + c \dfrac{\partial }{\partial x} \right) u(x,t) &= 0 \\ \end{align*} $$
このとき波動方程式の解$u(x,t)$は次の2つの方程式のいずれかを満たす。
$$ \left(\dfrac{\partial }{\partial t} - c \dfrac{\partial }{\partial x} \right) u(x,t) = 0 \tag{1} $$
$$ \left(\dfrac{\partial }{\partial t} + c \dfrac{\partial }{\partial x} \right) u(x,t) = 0 \tag{2} $$
$(1)$を満たす$u_{+}(x, t) = f(x + ct)$は$-x$方向に移動する波動である。反対に$(2)$を満たす$u_{-}(x, t) = g(x - ct)$は$+x$方向に移動する波動である。境界の$x = 0$の場所では$U$は左方向に行く波動だけが存在するはずなので、$U$は$(1)$を満たさなければならない。$x = 0$での$t$に対する微分係数を求めると、
$$ \dfrac{\partial u}{\partial t}|_{x = 0} \approx \dfrac{U_{0}^{n+1} - U_{0}^{n}}{\Delta t} \approx \dfrac{\dfrac{U_{0}^{n+1} + U_{1}^{n+1}}{2} - \dfrac{U_{0}^{n} + U_{1}^{n}}{2}}{\Delta t} $$
2つ目の近似式で関数値を隣接する2点の平均としたのは、より良い近似のためである2。同様に$x$に対する微分係数は次のように求めることができる。
$$ c\dfrac{\partial u}{\partial x}|_{x = 0} \approx c \dfrac{U_{1}^{n} - U_{0}^{n}}{\Delta x} \approx c\dfrac{\dfrac{U_{1}^{n+1} + U_{1}^{n}}{2} - \dfrac{U_{0}^{n+1} + U_{0}^{n}}{2}}{\Delta x} $$
$U$は$x = 0$で$(1)$を満たさなければならないので、
$$ \dfrac{\dfrac{U_{0}^{n+1} + U_{1}^{n+1}}{2} - \dfrac{U_{0}^{n} + U_{1}^{n}}{2}}{\Delta t} - c\dfrac{\dfrac{U_{1}^{n+1} + U_{1}^{n}}{2} - \dfrac{U_{0}^{n+1} + U_{0}^{n}}{2}}{\Delta x} = 0 $$
$$ \implies \dfrac{ U_{0}^{n+1} + U_{1}^{n+1} - U_{0}^{n} - U_{1}^{n}}{2 \Delta t} - c\dfrac{ U_{1}^{n+1} + U_{1}^{n} - U_{0}^{n+1} - U_{0}^{n}}{2 \Delta x} = 0 $$
$$ \implies U_{0}^{n+1} + U_{1}^{n+1} - U_{0}^{n} - U_{1}^{n} - \dfrac{c\Delta t}{\Delta x} \left( U_{1}^{n+1} + U_{1}^{n} - U_{0}^{n+1} - U_{0}^{n} \right) = 0 $$
$$ \implies \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{0}^{n+1} = \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{1}^{n} - (U_{1}^{n+1} - U_{0}^{n}) + \dfrac{c\Delta t}{\Delta x} \left( U_{1}^{n+1} - U_{0}^{n} \right) $$
$$ \implies \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{0}^{n+1} = \left(\dfrac{c \Delta t}{\Delta x} + 1 \right) U_{1}^{n} + \left(\dfrac{c \Delta t}{\Delta x} - 1 \right) \left( U_{1}^{n+1} - U_{0}^{n} \right) $$
$$ \begin{align*} \implies U_{0}^{n+1} &= U_{1}^{n} + \dfrac{\left(\dfrac{c \Delta t}{\Delta x} - 1 \right)}{\left(\dfrac{c \Delta t}{\Delta x} + 1 \right)} \left( U_{1}^{n+1} - U_{0}^{n} \right) \\ &= U_{1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{1}^{n+1} - U_{0}^{n} \right) \end{align*} $$
同様に、$x=1$の場合の吸収境界条件は次の通りである。
$$ U_{M}^{n+1} = U_{M-1}^{n} + \dfrac{\left(c \Delta t - \Delta x \right)}{\left(c \Delta t + \Delta x \right)} \left( U_{M-1}^{n+1} - U_{M}^{n} \right) $$
■
コード
1次元波動方程式における吸収境界条件をシミュレートしたコードは次の通りである。
using Plots
x = LinRange(-1, 1, 300)
Δx = x[2] - x[1]
t = LinRange(0, 4, 600)
Δt = t[2] - t[1]
μ = Δt / Δx
C = (μ - 1) / (μ + 1)
U = zeros(length(x), length(t))
U[:, 1] = exp.(-30 * x.^2)
U[:, 2] = exp.(-30 * x.^2)
for i ∈ 3:length(t) # without ABCs
U[2:end-1, i] = (μ^2)*U[3:end, i-1] + 2(1-μ^2)*U[2:end-1, i-1] + (μ^2)*U[1:end-2, i-1] - U[2:end-1, i-2]
end
V = zeros(length(x), length(t))
V[:, 1] = exp.(-30 * x.^2)
V[:, 2] = exp.(-30 * x.^2)
for i ∈ 3:length(t) # with ABCs
V[2:end-1, i] = (μ^2)*V[3:end, i-1] + 2(1-μ^2)*V[2:end-1, i-1] + (μ^2)*V[1:end-2, i-1] - V[2:end-1, i-2]
V[1, i] = V[2, i-1] + C*(V[2, i] - V[1, i-1]) # ABCs
V[end, i] = V[end-1, i-1] + C*(V[end-1, i] - V[end, i-1]) # ABCs
end
anim = @animate for i ∈ 1:2:Int(length(t)/3)
p1 = plot(x, U[:, i], xlims=(-1, 1), ylims=(-1,1), legend=false, dpi=300, title = "w/o ABCs")
p2 = plot(x, V[:, i], xlims=(-1, 1), ylims=(-1,1), legend=false, dpi=300, title = "w/ ABCs")
plot(p1, p2, layout=(2,1))
end
gif(anim, fps=40)
2次元での吸収境界条件をシミュレートしたコードは次の通りである。
using Plots
# 파라미터 설정
nx, ny = 200, 200 # 격자 크기
nt = 400 # 시간 스텝 수
Lx, Ly = 10.0, 10.0 # 공간 영역 크기
T = 10.0 # 총 시간
c = 1.0 # 파동 속도
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
dt = T / (nt - 1)
κ = c * dt / dx
# 초기 조건 설정
u = zeros(Float64, nx, ny, nt)
v = zeros(Float64, nx, ny, nt)
for i ∈ 1:nx, j ∈ 1:ny
if 13^2 < (i-160)^2 + (j-160)^2 < 15^2
u[i, j, 1] = 1
u[i, j, 2] = 1
v[i, j, 1] = 1
v[i, j, 2] = 1
end
if (i-40)^2 + (j-40)^2 < 4^2
u[i, j, 1] = 1
u[i, j, 2] = 1
v[i, j, 1] = 1
v[i, j, 2] = 1
end
end
for n in 2:(nt - 1)
u_xx = (u[3:end, 2:end-1, n] - 2 * u[2:end-1, 2:end-1, n] + u[1:end-2, 2:end-1, n]) / dx^2
u_yy = (u[2:end-1, 3:end, n] - 2 * u[2:end-1, 2:end-1, n] + u[2:end-1, 1:end-2, n]) / dy^2
u[2:end-1, 2:end-1, n + 1] = 2 * u[2:end-1, 2:end-1, n] - u[2:end-1, 2:end-1, n - 1] + c^2 * dt^2 * (u_xx + u_yy)
v_xx = (v[3:end, 2:end-1, n] - 2 * v[2:end-1, 2:end-1, n] + v[1:end-2, 2:end-1, n]) / dx^2
v_yy = (v[2:end-1, 3:end, n] - 2 * v[2:end-1, 2:end-1, n] + v[2:end-1, 1:end-2, n]) / dy^2
v[2:end-1, 2:end-1, n + 1] = 2 * v[2:end-1, 2:end-1, n] - v[2:end-1, 2:end-1, n - 1] + c^2 * dt^2 * (v_xx + v_yy)
u[1, 2:end-1, n+1] = u[2, 2:end-1, n] + ((κ-1)/(κ+1))*(u[2, 2:end-1, n] - u[1, 2:end-1, n-1])
u[end, 2:end-1, n+1] = u[end-1, 2:end-1, n] + ((κ-1)/(κ+1))*(u[end-1, 2:end-1, n] - u[end, 2:end-1, n-1])
u[2:end-1, 1, n+1] = u[2:end-1, 2, n] + ((κ-1)/(κ+1))*(u[2:end-1, 2, n] - u[2:end-1, 1, n-1])
u[2:end-1, end, n+1] = u[2:end-1, end-1, n] + ((κ-1)/(κ+1))*(u[2:end-1, end-1, n] - u[2:end-1, end, n-1])
end
b2r = cgrad([:blue, :white, :red])
# 애니메이션 생성
anim = @animate for n in 1:nt
p1 = heatmap(u[:, :, n], c=b2r, clim=(-1,1), legend=false, xticks=false, yticks=false, framestyle=:box, ratio=1, dpi=300, title="w/ ABCs")
p2 = heatmap(v[:, :, n], c=b2r, clim=(-1,1), legend=false, xticks=false, yticks=false, framestyle=:box, ratio=1, dpi=300, title="w/o ABCs")
plot(p1, p2, layout=(1,2), size=(728, 728/2))
end
gif(anim, fps=50)
Mur, Gerrit. “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations.” IEEE transactions on Electromagnetic Compatibility 4 (1981): 377-382. ↩︎