logo

Mur 흡수 경계 조건

Mur 흡수 경계 조건

Overview

Methods for numerically calculating the propagation of waves include the Finite Element Method (FEM) and the kk-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.

2t2u(x,t)=c22x2u(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)

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

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

Derivation1

Let’s reorganize the wave equation as follows.

2t2u(x,t)=c22x2u(x,t)    2t2u(x,t)c22x2u(x,t)=0    (tcx)(t+cx)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*}

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

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

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

(1)(1), satisfying u+(x,t)=f(x+ct)u_{+}(x, t) = f(x + ct), is a wave moving in the x-x direction. Conversely, u(x,t)=g(xct)u_{-}(x, t) = g(x - ct), satisfying (2)(2), is a wave moving in the +x+x direction. At the boundary x=0x = 0, only waves moving to the left should exist, so UU must satisfy (1)(1). Differentiating tt at x=0x = 0:

utx=0U0n+1U0nΔtU0n+1+U1n+12U0n+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}

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

cuxx=0cU1nU0nΔxcU1n+1+U1n2U0n+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}

Since UU must satisfy (1)(1) at x=0x = 0:

U0n+1+U1n+12U0n+U1n2ΔtcU1n+1+U1n2U0n+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+1U0nU1n2ΔtcU1n+1+U1nU0n+1U0n2Δ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+1U0nU1ncΔtΔx(U1n+1+U1nU0n+1U0n)=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+1U0n)+cΔtΔx(U1n+1U0n) \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Δx1)(U1n+1U0n) \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Δx1)(cΔtΔx+1)(U1n+1U0n)=U1n+(cΔtΔx)(cΔt+Δx)(U1n+1U0n) \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=1x=1 is as follows.

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

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