Mur 흡수 경계 조건
Overview
Methods for numerically calculating the propagation of waves include the Finite Element Method (FEM) and the -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.
Suppose this is simulated in a finite grid space . If we divide the -axis space into segments, it becomes , and if we divide time into segments, it becomes . The discretization of the solution of the wave equation is called . In this case, and . The absorbing boundary condition is as follows.
Explanation
The results of applying absorbing boundary conditions in 1D and 2D are as follows.
Derivation1
Let’s reorganize the wave equation as follows.
In this case, the solution of the wave equation satisfies one of the following two equations.
, satisfying , is a wave moving in the direction. Conversely, , satisfying , is a wave moving in the direction. At the boundary , only waves moving to the left should exist, so must satisfy . Differentiating at :
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 is obtained as follows.
Since must satisfy at :
Similarly, the absorbing boundary condition for is as follows.
■
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)
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. ↩︎