[PYTHON] Acoustic simulation FDTD method Mur absorption boundary derivation and implementation

Introduction

The FDTD method is one of the methods for numerical simulation of acoustics and electromagnetics. This is calculated by directly replacing the derivative of the PDE with the difference, which has the advantage of simplifying the theory and programming.

Since the FDTD method calculates a finite area, it is important to handle the boundary. If you leave the boundary at 0, it will be a fixed end and the wave will be totally reflected. If you don't want it to be reflected like an anechoic chamber, such as an open outdoor sound, you need to devise a boundary. This time, we will derive, implement, and try Mur's primary absorbing boundary, which is probably the simplest absorbing boundary.

I'm not specialized in this area, so it's quite possible that I'm wrong. At that time, I would appreciate it if you could let me know in the comments.

The original image doesn't have noise on the border like this, but when I use the online GIF compression tool, it looks like this.

reference

I referred to here for the discretization of the wave equation. Qiita: Numerical solution of wave equation

Refer to here for the absorbing boundary of Mur.

[Current state of electromagnetic field analysis method for light waves Time domain difference method-Expecting application to the optical field-](https://annex.jsap.or.jp/photonics/kogaku/public/27-11-kaisetsu5. pdf)

FDTD method as a photoelectric magnetic field analysis method

FDTD method

Discretization of wave equation

This time we will consider a two-dimensional wave equation. Where $ p $ is the pressure and $ c $ is the speed of sound.

\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2}- \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0

Discrete by center difference.

\frac{p^{n+1}_{i,j} - 2p^n_{i,j} + p^{n-1}_{i,j}}{\Delta t^2} = c^2 \left( \frac{p^{n}_{i+1,j} - 2p^n_{i,j} + p^{n}_{i-1,j}}{\Delta x^2} + \frac{p^{n}_{i,j+1} - 2p^n_{i,j} + p^{n}_{i,j-1}}{\Delta y^2} \right)

The formula is transformed so that $ p $ at the next time comes to the left side. Now let's assume that $ \ Delta x $ and $ \ Delta y $ are equal.

p^{n+1}_{i, j} = 2p^n_{i,j} - p^{n-1}_{i,j} + \frac{\Delta t^2 c^2}{\Delta x^2} \left( p^n_{i-1,j} + p^n_{i+1,j} + p^n_{i,j-1} + p^n_{i,j+1} - 4p^n_{i,j} \right)

FDTD program fixed end

The program was created and executed in Jupyter lab.

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

#The unit is m,s,m/s
#dx =Make dy(Because I did that with the formula transformation)
x_len = 8.0 #Length in x direction m
y_len = 8.0 #Length in y direction m
nx = 400 #Number of grids in the x direction
ny = 400 #Number of grids in the y direction
dx = x_len / nx #Grid width in x direction m
dy = y_len / ny #Lattice width m in the y direction
dt = 0.000001 #Time step s
c = 340.0 #Speed of sound m/s
C1 = (dt**2)*(c**2)/(dx**2)
C2 = (c*dt)/dx
x = np.linspace(0.0,x_len,nx)
y = np.linspace(0.0,y_len,ny)
X,Y = np.meshgrid(x,y)

if not dx == dy:
    print("Make dx and dy the same")

if c*dt/dx > 1.0:
    print("CFL={} < 1".format(c*dt/dx))

def get_nearest(xp,yp):
    """
For calculating the location of the sound source
    """
    return np.argmin((x-xp)**2),np.argmin((y-yp)**2)

sx,sy = get_nearest(4.0,4.0) #Sound source position

@jit
def sound_source(p,p_pre,t,f):
    p[sx,sy] = np.sin(2.0*np.pi*f*t)
    p_pre[sx,sy] = np.sin(2.0*np.pi*f*(t-dt))

@jit
def update(p,p_pre,t):
    #Vibration
    f = 1.0e3
    if t < 1.0/f * 5.0:
        sound_source(p,p_pre,t,f)

    p_next = np.zeros_like(p)
    for i in range(1,nx-1):
        for j in range(1,ny-1):
            p_next[i,j]  = 2.0*p[i,j] - p_pre[i,j] + C1*(p[i-1,j] + p[i+1,j] + p[i,j-1] + p[i,j+1] - 4.0*p[i,j])
    return p_next,p

p = np.zeros((nx,ny)) #Current time step
p_pre = np.zeros((nx,ny)) #Previous time step

t = 0
img_num = 1
for i in range(50000):
    if i%200==0:
        fig,axes = plt.subplots(figsize=(8,8))
        axes.imshow(np.flipud(p.T),vmin=-0.1,vmax=0.1,cmap="jet")
        fig.savefig("anime/{}.png ".format(img_num))
        plt.close()
        img_num += 1
    p,p_pre = update(p,p_pre,t)
    t += dt

I saved one graph in 0.2ms and animated it with GIMP.

Leaving the boundary at 0 results in fixed-end reflections.

Derivation of Mur's primary absorbing boundary

Consider a one-dimensional wave equation.

\frac{\partial^2 p}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0

Factor into the form $ a ^ 2-b ^ 2 = (a + b) (a-b) $.

\left( \frac{\partial p}{\partial x} + \frac{1}{c} \frac{\partial p}{\partial t} \right) \left( \frac{\partial p}{\partial x} - \frac{1}{c} \frac{\partial p}{\partial t} \right) = 0

The first term of the above equation seems to represent a receding wave, and it will be 0 if there is no reflection.

\frac{\partial p}{\partial x} + \frac{1}{c} \frac{\partial p}{\partial t} = 0

We will differentiate this. Space is at $ i- \ frac {1} {2} $ and time is at $ n + \ frac {1} {2} $.

\frac{p^{n+\frac{1}{2}}_{i} - p^{n+\frac{1}{2}}_{i-1}}{\Delta x} + \frac{1}{c} \frac{p^{n+1}_{i-\frac{1}{2}} - p^n_{i-\frac{1}{2}}}{\Delta t} = 0

Since there are no grid points at the positions of $ i- \ frac {1} {2} $ and $ n + \ frac {1} {2} $, the average value is used.

\frac{ \frac{p^{n+1}_{i}+p^n_{i}}{2} - \frac{p^{n+1}_{i-1}+p^n_{i-1}}{2}}{\Delta x} + \frac{1}{c} \frac{\frac{p^{n+1}_{i} + p^{n+1}_{i-1}}{2} - \frac{ p^n_{i} + p^n_{i-1}}{2}}{\Delta t} = 0

After that, transform the formula and bring $ p ^ {n + 1} _i $ to the left side and the others to the right side.

First, multiply both sides by $ 2c \ Delta t \ Delta x $

c \Delta t \left( p^{n+1}_i + p^n_i - p^{n+1}_{i-1} - p^{n+1}_{i-1} \right) + \Delta x \left( p^{n+1}_i + p^{n+1}_{i-1} - p^n_i - p^n_{i-1} \right) = 0

There are four types of $ p $, so I will summarize each one.

\left( c \Delta t + \Delta x \right) p^{n+1}_i
+ \left( c \Delta t - \Delta x \right) p^n_i
- \left( c \Delta t - \Delta x \right) p^{n+1}_{i-1}
- \left( c \Delta t + \Delta x \right) p^n_{i-1} = 0

Move all but the first term on the left side to the right side

\left( c \Delta t + \Delta x \right) p^{n+1}_i = 
 \left( c \Delta t - \Delta x \right) \left( p^{n+1}_{i-1} - p^n_i \right)
 + \left( c \Delta t + \Delta x \right) p^n_{i-1}

Divide both sides by $ \ left (c \ Delta t + \ Delta x \ right) $

 p^{n+1}_i = 
 \frac{\left( c \Delta t - \Delta x \right)}{\left( c \Delta t + \Delta x \right)} \left( p^{n+1}_{i-1} - p^n_i \right)
 + p^n_{i-1}

The formula has the same shape as the reference.

FDTD Program Mur Primary Absorption Boundary

Write only the changed part.

def update(p,p_pre,t):
    #Vibration
    f = 1.0e3
    if t < 1.0/f * 5.0:
        sound_source(p,p_pre,t,f)

    p_next = np.zeros_like(p)
    for i in range(1,nx-1):
        for j in range(1,ny-1):
            p_next[i,j]  = 2.0*p[i,j] - p_pre[i,j] + C1*(p[i-1,j] + p[i+1,j] + p[i,j-1] + p[i,j+1] - 4.0*p[i,j])

    #Mur's primary absorbing boundary
    for i in range(nx):
        p_next[i,0] = (c*dt-dx)/(c*dt+dx)*(p_next[i,1]-p[i,0]) + p[i,1]
        p_next[i,-1] = (c*dt-dx)/(c*dt+dx)*(p_next[i,-2]-p[i,-1]) + p[i,-2]

    for j in range(ny):
        p_next[0,j] = (c*dt-dx)/(c*dt+dx)*(p_next[1,j]-p[0,j]) + p[1,j]
        p_next[-1,j] = (c*dt-dx)/(c*dt+dx)*(p_next[-2,j]-p[-1,j]) + p[-2,j]

    return p_next,p

Waves that are vertically incident often disappear, but those that are obliquely incident remain a little. Boundary noise is due to GIF compression.

Derivation of Mur's secondary absorbing boundary

Mur's quadratic absorbing boundary can be derived by performing the same operation on the two-dimensional wave equation.

\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2}- \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2}=0

Since it is difficult to write a lot of partial differential symbols, the formula is transformed by setting the three terms as $ a $, $ b $, and $ c $, respectively.

Change the form of $ a ^ 2 + b ^ 2 --c ^ 2 = 0 $ to $ b ^ 2-\ left (c ^ 2 --a ^ 2 \ right) = 0 $

b^2 - \left[  c^2 \left\{ 1 -  \left( \frac{a}{c} \right) ^2 \right\} \right] = 0

\left( b + c \sqrt{1- \left( \frac{a}{c} \right) ^2} \right) \left( b - c \sqrt{1- \left( \frac{a}{c} \right) ^2} \right) = 0

The condition that there is no receding wave in the $ y $ direction

b + c \sqrt{1- \left( \frac{a}{c} \right) ^2} = 0

In the reference, the formula transformation ends at once from here, but it seems that the square root term is approximated by Taylor expansion.

When $ x \ ll 1 $, $ \ sqrt {1-x ^ 2} \ approx 1-\ frac {x ^ 2} {2} $

b + c \left\{ 1 - \frac{1}{2} \left( \frac{a}{c} \right) ^2\right\} = 0

Multiply $ c $ on both sides

bc + c^2 - \frac{1}{2} a^2 = 0

Put the sign back

\frac{\partial p}{\partial y \partial t} + \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} - \frac{c^2}{2} \frac{\partial^2 p}{\partial x^2} = 0

This is the conditional expression of Mur's secondary absorbing boundary.

When I tried to do the difference as well as the primary, it didn't work because the term of the future position such as $ p ^ {n + 1} _ {i-1, j} $ remained. Moreover, if you program the diff of the bibliography as it is, it will diverge and I am not sure.

Summary

One of the important absorbing boundaries was derived and implemented by the FDTD method. When dealing with phenomena that require more precision, such as diffraction, it is said to use a PML absorbing boundary with less reflection. If I have a chance, I would like to learn and post it.