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차원에서 흡수 경계 조건을 적용한 결과는 다음과 같다.

유도¹

파동 방정식을 다음과 같이 정리하자.

$$ \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)$는 다음의 두 방정식 중 하나를 만족한다.

$$ \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} $$

두번째 근사식에서 함숫값을 인접한 두 점의 평균으로 둔 것은 더 나은 근사를 위한 것이다². 마찬가지로 $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)