1 Acoustic wave equation

The acoustic wave equation in vector form is given by

$$ p_{tt}=c^2 \nabla^2 p \tag{1} $$

where, $p$ is the acoustic pressure, $c$ is the propagation speed of acoustic wave, $p_{tt}$ is the second time derivative of $p$, and $\nabla^2$ denotes the $Laplacian$ operator. We will discretize the equation using finite difference method in the following parts.

2 Propagation in 1D space

The Taylar’s series of $p(t+\Delta t)$ and $p(t-\Delta t)$ are given by

$$ \left \{ \begin{aligned} p(t+\Delta t) &= p(t) + p’(t)\Delta t + \frac{1}{2}p’'(t)\Delta t^2 + \mathcal{O}(\Delta t^2)\\
p(t-\Delta t) &= p(t) - p’(t)\Delta t + \frac{1}{2}p’'(t)\Delta t^2 + \mathcal{O}(\Delta t^2) \end{aligned} \right. \tag{2}, $$

so we get $p_{tt}$ in discrete form as

$$ p_{tt}=\frac{\partial ^2 p(x, t)}{\partial t^2}= \frac{p(x, t+\Delta t)-2p(x, t)+p(x, t-\Delta t)}{\Delta t^2} \tag{3}. $$

Similarly, the discretae form of $\nabla^2p$ becomes

$$ \nabla^2 p = \frac{p(x+\Delta x)-2p(x)+p(x-\Delta x)}{\Delta x^2} \tag{4}. $$

Lastly, we get the 3-point finite difference in time and space domains of 1D acoustic wave and the updated pressure is given by

$$ \begin{aligned} p(x, t+\Delta t) &= 2p(x, t) - p(x, t-\Delta t)\\
&+ \frac{c^2\Delta t^2}{\Delta x^2}[p(x+\Delta x, t)-2p(x, t)+p(x-\Delta x, t)] \end{aligned} \tag{5}. $$

For simplicity, we rewrite eq. (5) as $$ p^{i+1}_j = 2p^i_j-p^{i-1}_j + \frac{c^2\Delta t^2}{\Delta x^2}[p^i_{j+1}-2p^i_j+p^i_{j-1}] \tag{6}, $$

where, $i, j$ indicate the indices of time and space, respectively.

3 Propagation in 2D space

For the acoustic wave in 2D homogenous space domains, we have

$$ \frac{\partial ^2p(x, t)}{\partial t^2} = c^2 [\frac{\partial ^2p(x, t)}{\partial x^2} + \frac{\partial ^2p(x, t)}{\partial y^2}] \tag{7}. $$

Here, we use $i, j, k$ to indicate the indices of time, $x$ and $y$, respectively. Thus, $$ p^{i+1}_{j, k}=2p^i_{i, j}-p^{i-1}_{j,k}+c^2\Delta t^2 [\frac{p^i_{j, k+1}-2p^i_{j,k}+p^i_{j,k-1}}{\Delta x^2}+ \frac{p^i_{j+1, k}-2p^i_{j,k}+p^i_{j-1,k}}{\Delta y^2}+ ] \tag{8}. $$

4 Stability condition

$$ \frac{c \Delta t}{\Delta x} \le 1 \tag{9}. $$

5 Numerical examples

We use Ricker wavelet as the excitation of the acoustic wave equation, and the Ricker wavelet is given by $$ r(t)=\{1-2[\pi f_0 (t-t_0)]^2\} e^{-[\pi f_0 (t-t_0)]^2} \tag{10}, $$ where, $f_0$ is the center frequency, and $t_0$ is the shift time. The following figure shows the shape of Ricker wavelet waveform in time domain with given parameters. RICKER

5.1 1D

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

import numpy as np
import matplotlib.pyplot as plt

def aco_1d(nt, dt, nx, dx, sx, c, s):
	# Initilization of pressure matrix
    u = np.zeros((nt, nx))
    A = (c*dt/dx) ** 2
    for i in range(1, nt-1):
        for j in range(1, nx-1):
			# Set value of pressure as source excitation
            if j == sx:
                u[i+1, j] = s[i]
			# Update pressure
            else:
                u[i+1, j] = 2*u[i, j] - u[i-1, j] \
                          + A*(u[i, j+1]-2*u[i, j]+u[i, j-1])
    return u

def main():
	# No. of time, time inverval and acoustic wave speed
    nt = 1001; dt = 1e-3; c = 340
	# Center frequency and shift time of Ricker wavelet
    f0 = 20; t = np.arange(nt) * dt; t0 = nt//15*dt
    s = (1-2*(np.pi*f0*(t-t0))**2) * np.exp(-(np.pi*f0*(t-t0))**2)
    t -= t0
	# No. of space
    nx = 501; dx = 1
	# Position index of source
	sx = nx // 2
    x = np.arange(nx) * dx
    u = aco_1d(nt, dt, nx, dx, sx, c, s)
    plt.figure(figsize=(5, 7))
    ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=1, colspan=1)
    ax1.plot(t, s, lw=1, color='k')
    ax1.set_xlabel('Time (s)', fontsize=13)
    ax1.set_title('Ricker wavelet', fontsize=12)
    ax1.tick_params(labelsize=11)

    ax2 = plt.subplot2grid((3, 1), (1, 0), rowspan=2, colspan=1)
    for i in range(0, len(u[0]), 20):
        ax2.plot(t, u[:, i]*15+i, lw=1, color='#555555', zorder=0)
    ax2.scatter(0, sx*dx, marker='*', s=150,
                facecolor='r', zorder=1)
    ax2.set_xlabel('Time (s)', fontsize=13)
    ax2.set_ylabel('X (m)', fontsize=15) 
    ax2.set_xlabel('Time (s)', fontsize=15)
    ax2.tick_params(labelsize=11)
    plt.tight_layout()
    plt.show()

if __name__ == '__main__':
    main()

1D_HOMO

5.2 2D

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

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

#@jit(nopython=True)
def main():
    nx = 501
    ny = 501
    dx = 0.5
    nt = 401
    dt = 1e-3
    c = 340
    rx = 150
    ry = 150
    sx = 200
    sy = 250
    f0 = 45
    t0 = 0.02
    t = np.arange(nt) * dt
    sr = (1-2*(np.pi*f0*(t-t0))**2) * np.exp(-(np.pi*f0*(t-t0))**2)
    u1 = np.zeros((nx, ny))
    u2 = np.zeros((nx, ny))
    u3 = np.zeros((nx, ny))
    d = np.zeros((11, nt))
    A = (dt*c/dx) ** 2
    for i in range(1, nt-1):
        for j in range(1, ny-1):
            for k in range(1, nx-1):
                if j == sy and k == sx:
                    u3[sy, sx] = sr[i]
                else:
                    u3[j, k] = 2*u2[j, k] - u1[j, k] + \
                           A * (u2[j+1, k]+u2[j, k+1]+\
                           u2[j-1, k]+u2[j, k-1]-\
                           4*u2[j, k])
        u1 = u2.copy()
        u2 = u3.copy()
    return u1

if __name__ == '__main__':
    u1  = main()
    np.savetxt('wave.txt', u1)
	plt.figure(figsize=(5.5, 5))
	plt.pcolormesh(d, cmap='coolwarm')
	plt.scatter(200, 250, marker='*', s=100, facecolor='k')
	plt.xlabel('X (m)', fontsize=13)
	plt.ylabel('Y (m)', fontsize=13)
	plt.xticks(fontsize=13)
	plt.yticks(fontsize=13)
	plt.show()

2D-HOMO