Mur 흡수 경계 조건

Overview

Methods for numerically calculating the propagation of waves include the Finite Element Method (FEM) and the $k$-space method. These methods assume that waves propagate infinitely, but in actual simulations, waves propagate within a finite grid. This can cause reflections at the boundaries. To solve this, grids can be set much larger than the actual region of interest, or boundary conditions can be applied to prevent wave reflections. These boundary conditions are called Absorbing Boundary Conditions.

Formulation

Consider the following wave equation.

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

Suppose this is simulated in a finite grid space $[0, 1] \times [0, T]$. If we divide the $x$-axis space into $M$ segments, it becomes $\Delta x = 1/M$, and if we divide time into $N$ segments, it becomes $\Delta t = T/N$. The discretization of the solution $u$ of the wave equation is called $U_{i}^{n} = u(i \ast \Delta x, n \ast \Delta t)$. In this case, $i = 0, 1, \cdots, M$ and $n = 0, 1, \cdots, N$. The absorbing boundary condition is as follows.

$$ 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) $$

Explanation

The results of applying absorbing boundary conditions in 1D and 2D are as follows.

Derivation¹

Let’s reorganize the wave equation as follows.

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

In this case, the solution $u(x,t)$ of the wave equation satisfies one of the following two equations.

$$ \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)$, satisfying $u_{+}(x, t) = f(x + ct)$, is a wave moving in the $-x$ direction. Conversely, $u_{-}(x, t) = g(x - ct)$, satisfying $(2)$, is a wave moving in the $+x$ direction. At the boundary $x = 0$, only waves moving to the left should exist, so $U$ must satisfy $(1)$. Differentiating $t$ at $x = 0$:

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

In the second approximation, the function value is set to the average of the two adjacent points for a better approximation². Similarly, the derivative concerning $x$ is obtained as follows.

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

Since $U$ must satisfy $(1)$ at $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 $$

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

Similarly, the absorbing boundary condition for $x=1$ is as follows.

$$ 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) $$

■

Code

The code simulating the absorbing boundary condition for the 1D wave equation is as follows.

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)

The code simulating the absorbing boundary condition in 2D is as follows.

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)