logo

吸収境界条件

吸収境界条件

概要

波動の伝播を数値的に計算する方法には有限要素法(FDM)、kk-space methodなどがある。これらの方法は波動が無限に伝播することを仮定しているが、実際のシミュレーションでは有限の格子内で波動が伝播する。そのため、境界では波動が反射する問題が生じる。これを解決するためには、実際にシミュレーションしたい領域よりもはるかに広く格子を設定するか、境界条件をうまく設定して波動の反射が起こらないようにする。そのような境界条件を吸収境界条件absorbing boundary condition (ABC)という。

公式

次のような波動方程式を考えよう。

∂2∂t2u(x,t)=c2∂2∂x2u(x,t)for (x,t)∈R×[0,∞) \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]×[0,T][0, 1] \times [0, T]でシミュレーションするとしよう。xx軸空間をMM等分するとΔx=1/M\Delta x = 1/M、時間をNN等分するとΔt=T/N\Delta t = T/Nである。波動方程式の解uuの離散化をUin=u(i∗Δx,n∗Δt)U_{i}^{n} = u(i \ast \Delta x, n \ast \Delta t)という。このときi=0,1,⋯ ,Mi = 0, 1, \cdots, M、n=0,1,⋯ ,Nn = 0, 1, \cdots, Nである。吸収境界条件は次のようになる。

U0n+1=U1n+(cΔt−Δx)(cΔt+Δx)(U1n+1−U0n) 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)

UMn+1=UM−1n+(cΔt−Δx)(cΔt+Δx)(UM−1n+1−UMn) 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

波動方程式を次のように整理しよう。

∂2∂t2u(x,t)=c2∂2∂x2u(x,t)  ⟹  ∂2∂t2u(x,t)−c2∂2∂x2u(x,t)=0  ⟹  (∂∂t−c∂∂x)(∂∂t+c∂∂x)u(x,t)=0 \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)u(x,t)は次の2つの方程式のいずれかを満たす。

(∂∂t−c∂∂x)u(x,t)=0(1) \left(\dfrac{\partial }{\partial t} - c \dfrac{\partial }{\partial x} \right) u(x,t) = 0 \tag{1}

(∂∂t+c∂∂x)u(x,t)=0(2) \left(\dfrac{\partial }{\partial t} + c \dfrac{\partial }{\partial x} \right) u(x,t) = 0 \tag{2}

(1)(1)を満たすu+(x,t)=f(x+ct)u_{+}(x, t) = f(x + ct)は−x-x方向に移動する波動である。反対に(2)(2)を満たすu−(x,t)=g(x−ct)u_{-}(x, t) = g(x - ct)は+x+x方向に移動する波動である。境界のx=0x = 0の場所ではUUは左方向に行く波動だけが存在するはずなので、UUは(1)(1)を満たさなければならない。x=0x = 0でのttに対する微分係数を求めると、

∂u∂t∣x=0≈U0n+1−U0nΔt≈U0n+1+U1n+12−U0n+U1n2Δ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。同様にxxに対する微分係数は次のように求めることができる。

c∂u∂x∣x=0≈cU1n−U0nΔx≈cU1n+1+U1n2−U0n+1+U0n2Δ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}

UUはx=0x = 0で(1)(1)を満たさなければならないので、

U0n+1+U1n+12−U0n+U1n2Δt−cU1n+1+U1n2−U0n+1+U0n2Δx=0 \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

  ⟹  U0n+1+U1n+1−U0n−U1n2Δt−cU1n+1+U1n−U0n+1−U0n2Δ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

  ⟹  U0n+1+U1n+1−U0n−U1n−cΔtΔx(U1n+1+U1n−U0n+1−U0n)=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

  ⟹  (cΔtΔx+1)U0n+1=(cΔtΔx+1)U1n−(U1n+1−U0n)+cΔtΔx(U1n+1−U0n) \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)

  ⟹  (cΔtΔx+1)U0n+1=(cΔtΔx+1)U1n+(cΔtΔx−1)(U1n+1−U0n) \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)

  ⟹  U0n+1=U1n+(cΔtΔx−1)(cΔtΔx+1)(U1n+1−U0n)=U1n+(cΔt−Δx)(cΔt+Δx)(U1n+1−U0n) \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=1x=1の場合の吸収境界条件は次の通りである。

UMn+1=UM−1n+(cΔt−Δx)(cΔt+Δx)(UM−1n+1−UMn) 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)

  1. 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. ↩︎

  2. https://eecs.wsu.edu/~schneidj/ufdtd/chap6.pdf ↩︎