# An Unconditionally Stable Method for Solving the Acoustic Wave Equation

Authors: Zhi-Kai Fu; Li-Hua Shi; Zheng-Yu Huang; Shang-Chen Fu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

An unconditionally stable method for solving the time-domain acoustic wave equation using Associated Hermit orthogonal functions is proposed. The second-order time derivatives in acoustic wave equation are expanded by these orthogonal basis functions. By applying Galerkin temporal testing procedure, the time variable can be eliminated from the calculations. The restriction of Courant-Friedrichs-Levy (CFL) condition in selecting time step for analyzing thin layer can be avoided. Numerical results show the accuracy and the efficiency of the proposed method.


---

## Body

## 1. Introduction

Numerical simulation of the acoustic wave equation has been widely used in many areas, such as geophysics, exploration, and ultrasonic detection [1–5]. The finite difference (FD) method in time domain is one of the most widely used methods in solving problems related to acoustics due to its efficiency and simplicity [6–8]. However, as the FD method is an explicit time-marching algorithm, the time step must be limited by the Courant-Friedrichs-Levy (CFL) condition, which means that the computation time will increase dramatically as the spatial grid becomes small. Hence, some unconditionally stable methods were proposed to solve this problem [9–12]. The main categories of these methods include the alternating-direction implicit (ADI) method [10] and the orthogonal decomposition method using weighted Laguerre polynomials (WLP) [11]. Recently, Huang et al. [12, 13] developed a new orthogonal decomposition method using Associated Hermit orthogonal functions and have applied it to first-order Maxwell’s equation.In this paper, we extended the unconditionally stable method using AH functions to solve the second-order acoustic wave equation. Firstly, the second-order time derivatives in the acoustic wave equation are expanded by a set of orthogonal basis functions. Secondly, the time variable can be eliminated from the calculations by using Galerkin’s method. Finally, a set of implicit equations in the whole computational domain is established and solved in the AH domain.


## 2. Numerical Formulation for the Acoustic Wave Equation

### 2.1. 1D Acoustic Wave Equation

For simplicity, we consider the 1D acoustic wave equation with an external source in isotropic media(1)
∂
2
U
∂
t
2
=
v
2
∂
2
U
∂
x
2
+
f
x
,
t
,
0
<
x
<
l
,
0
<
t
<
T
,
U
x
,
0
=
0
,
∂
U
x
,
0
∂
t
=
0
,where v is the sound velocity, U is the wave displacement, and f is an external source. In order to truncate the computational domain, we choose Mur’s first-order absorbing boundary condition (ABC) for infinite problem from [14]. At the boundary points x
0
=
0 or x
L
=
l, we have(2)
∂
U
∂
x
-
1
c
∂
U
∂
t
x
=
0
or
l
=
0
.Consider a set of modified AH orthogonal basis functions given by [12](3)
ϕ
n
t
=
2
n
n
!
σ
π
1
/
2
-
1
/
2
e
-
t
-
2
/
2
-
1
n
e
t
-
2
d
n
/
d
t
n
e
-
t
-
2
,where t
¯
=
(
t
-
T
f
)
/
σ is the transformed time variable, T
f is the time-translating parameter, and σ is the time-scaling parameter. By selecting proper T
f and σ, the modified AH orthogonal basis functions can be used to approximate the causal displacement U
(
x
,
t
) by(4)
U
x
,
t
=
∑
n
=
0
N
-
1
u
n
ϕ
n
t
,where u
n
=
∫
-
∞
∞
U
(
x
,
t
)
φ
n
(
t
)
d
t is the nth expansion coefficient; N is selected according to the bandwidth of the analyzed signal [15].With the time derivation property of AH functions [15](5)
d
d
t
φ
n
t
=
-
1
σ
1
2
φ
1
t
n
=
0
1
σ
n
2
φ
n
-
1
t
-
1
σ
n
+
1
2
φ
n
+
1
t
n
≥
1
.We can express the first-order time derivation and the second-order time derivation of the displacementU
(
x
,
t
) as(6)
∂
U
x
,
t
∂
t
=
∑
n
=
0
∞
u
n
+
1
x
n
+
1
2
σ
2
-
u
n
-
1
x
n
2
σ
2
·
ϕ
n
t
,
(7)
∂
2
∂
t
2
U
x
,
t
=
∑
n
=
0
∞
u
n
+
2
x
n
+
2
2
σ
2
n
+
1
2
σ
2
-
2
n
+
1
2
σ
2
u
n
x
+
u
n
-
2
x
n
-
1
2
σ
2
n
2
σ
2
ϕ
n
t
,where u
-
2
(
x
)
=
0, u
-
1
(
x
)
=
0, and n represents the order in the AH domain. The deduction can be found in the Appendix.Substituting (4) and (7) into (1), we obtain(8)
∑
n
=
0
∞
u
n
+
2
x
n
+
2
2
σ
2
n
+
1
2
σ
2
-
2
n
+
1
2
σ
2
u
n
x
+
u
n
-
2
x
n
-
1
2
σ
2
n
2
σ
2
ϕ
n
t
=
v
2
∑
n
=
0
∞
∂
2
∂
x
2
u
n
x
·
ϕ
n
t
+
∑
n
=
0
∞
f
n
x
ϕ
n
t
.Based on the orthogonal property of the AH functions, we can obtain equations of the AH expansion coefficients using the temporal Galerkin testing procedure [11](9)
u
n
+
2
x
n
+
2
2
n
+
1
2
-
2
n
+
1
2
u
n
x
+
u
n
-
2
x
n
-
1
2
n
2
=
v
2
σ
2
∂
2
∂
x
2
u
n
x
+
σ
2
f
n
x
,where f
n
(
x
)
=
∫
-
T
s
/
2
T
s
/
2
f
(
x
,
t
)
φ
n
(
t
)
d
t, n
=
0,1
,
2
,
…. We can partition the spatial domain [
0
,
l
] into 0
=
x
1
<
x
1
<
·
·
·
<
x
L
=
l. Then the spatial discrete form of (9) at point x
i
(
0
<
i
<
L
) can be written as(10)
u
i
n
+
2
n
+
2
2
n
+
1
2
-
2
n
+
1
2
u
i
n
+
u
i
n
-
2
n
-
1
2
n
2
=
v
2
σ
2
Δ
x
2
u
i
+
1
n
-
2
u
i
n
+
u
i
-
1
n
+
σ
2
f
i
n
,where Δ
x is the spatial grid size.In (10), a set of implicit equations is included in the AH domain. Based on the recurrence relation from the n
-
2th order coefficients and the nth order coefficients to the n
+
2th order coefficients, we can rewrite (10) in a matrix form:(11)
-
v
2
σ
2
Δ
x
2
β
U
i
+
1
+
α
U
i
-
v
2
σ
2
Δ
x
2
β
U
i
-
1
=
σ
2
β
f
i
,where [
β
]
=
I
N
×
N in a unit matrix. Consider(12)
α
=
-
1
2
-
2
v
2
σ
2
Δ
x
2
2
2
1
2
-
3
2
-
2
v
2
σ
2
Δ
x
2
3
2
2
2
1
2
2
2
-
5
2
-
2
v
2
σ
2
Δ
x
2
⋱
⋱
⋱
N
-
2
2
N
-
1
2
-
2
N
-
1
+
1
2
-
2
v
2
σ
2
Δ
x
2
N
×
N
.Here,[
α
] is a tridiagonal coefficient matrix. [
U
]
i
=
(
u
i
0
,
u
i
1
,
…
,
u
i
N
-
1
)
T and [
f
]
i
=
(
f
i
0
,
f
i
1
,
…
,
f
i
N
-
1
)
T are the N-tuple representations in AH subspace for the displacement and the external source, respectively. In (11), we can find each point has a relationship with its adjacent two points, and a set of implicit equations can be obtained in the whole computation domain from point x
1 to x
L:(13)
A
U
=
J
,where(14)
A
=
α
-
v
2
σ
2
Δ
x
2
β
-
v
2
σ
2
Δ
x
2
β
α
-
v
2
σ
2
Δ
x
2
β
-
v
2
σ
2
Δ
x
2
β
α
-
v
2
σ
2
Δ
x
2
β
⋱
⋱
⋱
-
v
2
σ
2
Δ
x
2
β
α
L
×
Lis a tridiagonal matrix with submatrix elements. [
U
] is a combination of all [
U
]
i in the whole computation domain. [
J
] is a vector related with external source expansion coefficients for all points.Substitute (4) and (6) into (2) and eliminate the temporal terms to obtain (at x
=
x
L)(15)
∂
∂
x
U
n
x
-
1
c
u
n
+
1
x
n
+
1
2
σ
2
-
u
n
-
1
x
n
2
σ
2
=
0
.Using the average techniques, we have(16)
U
n
L
-
1
2
=
U
n
L
-
1
+
U
n
L
2
,
∂
∂
x
U
n
L
-
1
2
=
U
n
L
-
U
n
L
-
1
Δ
x
.Substituting (16) into (15), we have(17)
2
c
σ
Δ
x
U
L
n
-
n
+
1
2
U
L
n
+
1
+
n
2
U
L
n
-
1
=
2
c
σ
Δ
x
U
L
-
1
n
+
n
+
1
2
U
L
-
1
n
+
1
-
n
2
U
L
-
1
n
-
1
.Applying (17) into all orders in the AH domain, we can obtain(18)
2
v
σ
Δ
x
β
+
ε
U
L
-
2
v
σ
Δ
x
β
-
ε
U
L
-
1
=
0
,where(19)
ε
=
-
1
2
1
2
⋱
⋱
-
N
-
1
2
N
-
1
2
N
×
N
.Using the technique above, the absorbing boundary condition at pointx
=
x
1 can be expressed as(20)
2
v
σ
Δ
x
β
+
ε
U
1
-
2
v
σ
Δ
x
β
-
ε
U
2
=
0
.Introducing the boundary condition (18) and (20) into (11), (13) can be replaced by(21)
A
-
U
=
J
-
,where [
A
-
] is an adaptive coefficient matrix including absorbing boundary condition and [
J
-
] is a coefficient matrix corresponding to [
A
-
].We can use the lower-upper (LU) decomposition method to decompose[
A
-
] and use the back-substitution method to calculate the coefficients of the displacement for all points. The displacement U
(
x
,
t
) can be reconstructed finally by using (4).


### 2.2. 2D Acoustic Wave Equation

Considering the need for practical application, we extend the method to the 2D acoustic wave equation, which is given by(22)
∂
2
U
∂
t
2
=
v
2
∂
2
U
∂
x
2
+
∂
2
U
∂
y
2
+
f
x
,
y
,
t
,
x
,
y
∈
Ω
×
0
,
T
,
U
x
,
y
,
0
=
0
,
∂
U
x
,
y
,
0
∂
t
=
0
,where Ω
=
[
0
,
l
1
]
×
[
0
,
l
2
] is divided into N
x
×
N
y grid.Substituting (4) and (7) into (22) and rearranging it, we can obtain a matrix form equation which is similar to (11):(23)
-
v
2
σ
2
Δ
x
2
β
U
i
-
1
,
j
-
v
2
σ
2
Δ
x
2
β
U
i
+
1
,
j
+
α
U
i
,
j
-
v
2
σ
2
Δ
y
2
β
U
i
,
j
-
1
-
v
2
σ
2
Δ
y
2
β
U
i
,
j
-
1
=
v
2
σ
2
β
f
i
,
j
,where Δ
x is the spatial grid size in x direction and Δ
y is the spatial grid size in y direction. In (23), we can find each point has a relationship with its adjacent four points, so matrix [
A
] is a five-diagonal coefficient matrix in 2D case.For the absorbing boundary condition, the boundary condition can be obtained by using (15) to (20). But the corner-point needs special treatment and the spatial grid size is adjusted to h
=
Δ
x
2
+
Δ
y
2 at corner-point. For the four corner-points, Mur’s absorbing boundary condition can be expressed as(24)
2
v
σ
h
β
+
ε
U
1.1
-
2
v
σ
h
β
-
ε
U
2.2
=
0
,
2
v
σ
h
β
+
ε
U
1
,
N
y
-
2
v
σ
h
β
-
ε
U
2
,
N
y
-
1
=
0
,
2
v
σ
h
β
+
ε
U
N
x
,
1
-
2
v
σ
h
β
-
ε
U
N
x
-
1,2
=
0
,
2
v
σ
h
β
+
ε
U
N
x
,
N
y
-
2
v
σ
h
β
-
ε
U
N
x
-
1
,
N
y
-
1
=
0
.Introducing Mur’s absorbing boundary condition, the coefficient matrix[
A
] can be adjusted to [
A
-
]. Using the lower-upper (LU) decomposition method and the back-substitution method, the coefficients of all points can be calculated. And the displacement U
(
x
,
y
,
t
) can be reconstructed finally using the method described in 1D case.


## 2.1. 1D Acoustic Wave Equation

For simplicity, we consider the 1D acoustic wave equation with an external source in isotropic media(1)
∂
2
U
∂
t
2
=
v
2
∂
2
U
∂
x
2
+
f
x
,
t
,
0
<
x
<
l
,
0
<
t
<
T
,
U
x
,
0
=
0
,
∂
U
x
,
0
∂
t
=
0
,where v is the sound velocity, U is the wave displacement, and f is an external source. In order to truncate the computational domain, we choose Mur’s first-order absorbing boundary condition (ABC) for infinite problem from [14]. At the boundary points x
0
=
0 or x
L
=
l, we have(2)
∂
U
∂
x
-
1
c
∂
U
∂
t
x
=
0
or
l
=
0
.Consider a set of modified AH orthogonal basis functions given by [12](3)
ϕ
n
t
=
2
n
n
!
σ
π
1
/
2
-
1
/
2
e
-
t
-
2
/
2
-
1
n
e
t
-
2
d
n
/
d
t
n
e
-
t
-
2
,where t
¯
=
(
t
-
T
f
)
/
σ is the transformed time variable, T
f is the time-translating parameter, and σ is the time-scaling parameter. By selecting proper T
f and σ, the modified AH orthogonal basis functions can be used to approximate the causal displacement U
(
x
,
t
) by(4)
U
x
,
t
=
∑
n
=
0
N
-
1
u
n
ϕ
n
t
,where u
n
=
∫
-
∞
∞
U
(
x
,
t
)
φ
n
(
t
)
d
t is the nth expansion coefficient; N is selected according to the bandwidth of the analyzed signal [15].With the time derivation property of AH functions [15](5)
d
d
t
φ
n
t
=
-
1
σ
1
2
φ
1
t
n
=
0
1
σ
n
2
φ
n
-
1
t
-
1
σ
n
+
1
2
φ
n
+
1
t
n
≥
1
.We can express the first-order time derivation and the second-order time derivation of the displacementU
(
x
,
t
) as(6)
∂
U
x
,
t
∂
t
=
∑
n
=
0
∞
u
n
+
1
x
n
+
1
2
σ
2
-
u
n
-
1
x
n
2
σ
2
·
ϕ
n
t
,
(7)
∂
2
∂
t
2
U
x
,
t
=
∑
n
=
0
∞
u
n
+
2
x
n
+
2
2
σ
2
n
+
1
2
σ
2
-
2
n
+
1
2
σ
2
u
n
x
+
u
n
-
2
x
n
-
1
2
σ
2
n
2
σ
2
ϕ
n
t
,where u
-
2
(
x
)
=
0, u
-
1
(
x
)
=
0, and n represents the order in the AH domain. The deduction can be found in the Appendix.Substituting (4) and (7) into (1), we obtain(8)
∑
n
=
0
∞
u
n
+
2
x
n
+
2
2
σ
2
n
+
1
2
σ
2
-
2
n
+
1
2
σ
2
u
n
x
+
u
n
-
2
x
n
-
1
2
σ
2
n
2
σ
2
ϕ
n
t
=
v
2
∑
n
=
0
∞
∂
2
∂
x
2
u
n
x
·
ϕ
n
t
+
∑
n
=
0
∞
f
n
x
ϕ
n
t
.Based on the orthogonal property of the AH functions, we can obtain equations of the AH expansion coefficients using the temporal Galerkin testing procedure [11](9)
u
n
+
2
x
n
+
2
2
n
+
1
2
-
2
n
+
1
2
u
n
x
+
u
n
-
2
x
n
-
1
2
n
2
=
v
2
σ
2
∂
2
∂
x
2
u
n
x
+
σ
2
f
n
x
,where f
n
(
x
)
=
∫
-
T
s
/
2
T
s
/
2
f
(
x
,
t
)
φ
n
(
t
)
d
t, n
=
0,1
,
2
,
…. We can partition the spatial domain [
0
,
l
] into 0
=
x
1
<
x
1
<
·
·
·
<
x
L
=
l. Then the spatial discrete form of (9) at point x
i
(
0
<
i
<
L
) can be written as(10)
u
i
n
+
2
n
+
2
2
n
+
1
2
-
2
n
+
1
2
u
i
n
+
u
i
n
-
2
n
-
1
2
n
2
=
v
2
σ
2
Δ
x
2
u
i
+
1
n
-
2
u
i
n
+
u
i
-
1
n
+
σ
2
f
i
n
,where Δ
x is the spatial grid size.In (10), a set of implicit equations is included in the AH domain. Based on the recurrence relation from the n
-
2th order coefficients and the nth order coefficients to the n
+
2th order coefficients, we can rewrite (10) in a matrix form:(11)
-
v
2
σ
2
Δ
x
2
β
U
i
+
1
+
α
U
i
-
v
2
σ
2
Δ
x
2
β
U
i
-
1
=
σ
2
β
f
i
,where [
β
]
=
I
N
×
N in a unit matrix. Consider(12)
α
=
-
1
2
-
2
v
2
σ
2
Δ
x
2
2
2
1
2
-
3
2
-
2
v
2
σ
2
Δ
x
2
3
2
2
2
1
2
2
2
-
5
2
-
2
v
2
σ
2
Δ
x
2
⋱
⋱
⋱
N
-
2
2
N
-
1
2
-
2
N
-
1
+
1
2
-
2
v
2
σ
2
Δ
x
2
N
×
N
.Here,[
α
] is a tridiagonal coefficient matrix. [
U
]
i
=
(
u
i
0
,
u
i
1
,
…
,
u
i
N
-
1
)
T and [
f
]
i
=
(
f
i
0
,
f
i
1
,
…
,
f
i
N
-
1
)
T are the N-tuple representations in AH subspace for the displacement and the external source, respectively. In (11), we can find each point has a relationship with its adjacent two points, and a set of implicit equations can be obtained in the whole computation domain from point x
1 to x
L:(13)
A
U
=
J
,where(14)
A
=
α
-
v
2
σ
2
Δ
x
2
β
-
v
2
σ
2
Δ
x
2
β
α
-
v
2
σ
2
Δ
x
2
β
-
v
2
σ
2
Δ
x
2
β
α
-
v
2
σ
2
Δ
x
2
β
⋱
⋱
⋱
-
v
2
σ
2
Δ
x
2
β
α
L
×
Lis a tridiagonal matrix with submatrix elements. [
U
] is a combination of all [
U
]
i in the whole computation domain. [
J
] is a vector related with external source expansion coefficients for all points.Substitute (4) and (6) into (2) and eliminate the temporal terms to obtain (at x
=
x
L)(15)
∂
∂
x
U
n
x
-
1
c
u
n
+
1
x
n
+
1
2
σ
2
-
u
n
-
1
x
n
2
σ
2
=
0
.Using the average techniques, we have(16)
U
n
L
-
1
2
=
U
n
L
-
1
+
U
n
L
2
,
∂
∂
x
U
n
L
-
1
2
=
U
n
L
-
U
n
L
-
1
Δ
x
.Substituting (16) into (15), we have(17)
2
c
σ
Δ
x
U
L
n
-
n
+
1
2
U
L
n
+
1
+
n
2
U
L
n
-
1
=
2
c
σ
Δ
x
U
L
-
1
n
+
n
+
1
2
U
L
-
1
n
+
1
-
n
2
U
L
-
1
n
-
1
.Applying (17) into all orders in the AH domain, we can obtain(18)
2
v
σ
Δ
x
β
+
ε
U
L
-
2
v
σ
Δ
x
β
-
ε
U
L
-
1
=
0
,where(19)
ε
=
-
1
2
1
2
⋱
⋱
-
N
-
1
2
N
-
1
2
N
×
N
.Using the technique above, the absorbing boundary condition at pointx
=
x
1 can be expressed as(20)
2
v
σ
Δ
x
β
+
ε
U
1
-
2
v
σ
Δ
x
β
-
ε
U
2
=
0
.Introducing the boundary condition (18) and (20) into (11), (13) can be replaced by(21)
A
-
U
=
J
-
,where [
A
-
] is an adaptive coefficient matrix including absorbing boundary condition and [
J
-
] is a coefficient matrix corresponding to [
A
-
].We can use the lower-upper (LU) decomposition method to decompose[
A
-
] and use the back-substitution method to calculate the coefficients of the displacement for all points. The displacement U
(
x
,
t
) can be reconstructed finally by using (4).


## 2.2. 2D Acoustic Wave Equation

Considering the need for practical application, we extend the method to the 2D acoustic wave equation, which is given by(22)
∂
2
U
∂
t
2
=
v
2
∂
2
U
∂
x
2
+
∂
2
U
∂
y
2
+
f
x
,
y
,
t
,
x
,
y
∈
Ω
×
0
,
T
,
U
x
,
y
,
0
=
0
,
∂
U
x
,
y
,
0
∂
t
=
0
,where Ω
=
[
0
,
l
1
]
×
[
0
,
l
2
] is divided into N
x
×
N
y grid.Substituting (4) and (7) into (22) and rearranging it, we can obtain a matrix form equation which is similar to (11):(23)
-
v
2
σ
2
Δ
x
2
β
U
i
-
1
,
j
-
v
2
σ
2
Δ
x
2
β
U
i
+
1
,
j
+
α
U
i
,
j
-
v
2
σ
2
Δ
y
2
β
U
i
,
j
-
1
-
v
2
σ
2
Δ
y
2
β
U
i
,
j
-
1
=
v
2
σ
2
β
f
i
,
j
,where Δ
x is the spatial grid size in x direction and Δ
y is the spatial grid size in y direction. In (23), we can find each point has a relationship with its adjacent four points, so matrix [
A
] is a five-diagonal coefficient matrix in 2D case.For the absorbing boundary condition, the boundary condition can be obtained by using (15) to (20). But the corner-point needs special treatment and the spatial grid size is adjusted to h
=
Δ
x
2
+
Δ
y
2 at corner-point. For the four corner-points, Mur’s absorbing boundary condition can be expressed as(24)
2
v
σ
h
β
+
ε
U
1.1
-
2
v
σ
h
β
-
ε
U
2.2
=
0
,
2
v
σ
h
β
+
ε
U
1
,
N
y
-
2
v
σ
h
β
-
ε
U
2
,
N
y
-
1
=
0
,
2
v
σ
h
β
+
ε
U
N
x
,
1
-
2
v
σ
h
β
-
ε
U
N
x
-
1,2
=
0
,
2
v
σ
h
β
+
ε
U
N
x
,
N
y
-
2
v
σ
h
β
-
ε
U
N
x
-
1
,
N
y
-
1
=
0
.Introducing Mur’s absorbing boundary condition, the coefficient matrix[
A
] can be adjusted to [
A
-
]. Using the lower-upper (LU) decomposition method and the back-substitution method, the coefficients of all points can be calculated. And the displacement U
(
x
,
y
,
t
) can be reconstructed finally using the method described in 1D case.


## 3. Numerical Results

### 3.1. 1D Case

In order to validate the proposed method, a 1D structure constituted by two isotropic media separated by a thin layer (d
=
0.1 m, v
=
5.0 km/s) is analyzed. As shown in Figure 1, the thin layer is located at 60 m, and the total simulation length is 104 m. Displacement U
(
x
,
t
) at two points p
1 (38 m) and p
2 (78.1 m) is observed. The sound velocity in media A is 2.5 km/s, and the sound velocity in media B is 3.7 km/s. Nonuniform spatial grids are used. In the thin layer Δ
x
=
0.01 m; in other parts Δ
x
=
1 m.Figure 1
Configuration of a thin layer (d
=
0.1 m) embedded in two isotropic media.The external source is located atx
=
26 m and its expression is f
(
t
)
=
-
5.76
f
0
2
[
1
-
16
0.6
f
0
t
-
1
)
2
exp
⁡
[
-
8
(
0.6
f
0
t
-
1
)
2
], where f
0
=
100 Hz. Time support for the simulation is 0.16 s. To expand f
(
t
) in AH domain, we select N
=
52, σ
=
0.008, and T
f
=
0.5. Due to the limitation of the CFL condition, the time step is set as 2
×
1
0
-
6 s for the FD method under the condition of Δ
x
=
0.01 m. Time step is not involved in AH domain computation. However, for AH decomposition, we still need to specify a sample interval; the time step Δ
t
=
0.002 s is good enough. The absolute error R is defined as R
=
U
a
(
t
)
-
U
c
(
t
) for comparison of the two methods. U
c
(
t
) is the reference data obtained from the FD method with Δ
x
=
0.01 m, and U
a
(
t
) is the numerical result obtained from the FD method (Δ
x
=
0.1 m) and the proposed method (Δ
x
=
0.01 m). The numerical results at p
1 and their errors are shown in Figure 2(a). Figure 2(b) is the numerical results at p
2. From Figure 2, it can be seen that the proposed method has a good agreement with the FD method. The absolute error of the proposed method (Δ
x
=
0.01 m) is much lower than that of the FD method (Δ
x
=
0.1 m) at both p
1 and p
2.Figure 2
The numerical results and their errors at (a)p
1 and (b) p
2.
(a)
(b)Table1 shows the comparison of CPU time between the proposed method and the FD method. Under the condition of Δ
x
=
0.01 m, the time step of the proposed method is 1000 times as large as that of the FD method, and the total CPU time for the proposed method can be reduced to about 10% of the FD method.Table 1
The comparison of CPU time.
Δ
x = 0.01 m
Δ
t
CPU time
FD method
2
×
1
0
-
6 s
3.198 s
Proposed method
2
×
1
0
-
3 s
0.3127 s


### 3.2. 2D Case

In this section, we show results from a 2D numerical modeling by the FD method and the proposed method. The 2D model is established by two isotropic media separated by a thin layer (d
=
0.01 m, v
=
6.0 km/s), which is similar to the 1D case. As shown in Figure 3, the computational domain is 60 m × 90 m, and the thin layer is located at 72 m in x direction. Displacement U
(
x
,
y
,
t
) at two points p
1
(
36
m
,
30
m
) and p
2
(
81.01
m
,
30
m
) is observed. The sound velocity in media A is 3 km/s, and the sound velocity in media B is 4.7 km/s. Nonuniform spatial grids are used. In x direction, Δ
x
=
0.001 m in the thin layer, and Δ
x
=
3 m in other parts. In y direction, Δ
y
=
3 m. The external source is the same as the 1D case and is located at s
(
27
m
,
30
m
).Figure 3
Configuration of a thin layer (d
=
0.01 m) embedded in two isotropic media.Time support for the simulation is 0.06 s. The time step is set as1.33
×
1
0
-
7 s and 5
×
1
0
-
4 s for the FD method and the proposed method, respectively. Figure 4 shows the 2D modeling records at p
1 and p
2 for the multilayer isotropic media model. The waveform corresponding to the proposed method is in good agreement with that by the FD method, which demonstrates that the proposed method has good accuracy. Table 2 shows the comparison of CPU time between the proposed method and the FD method. Under the condition of Δ
x
=
0.001 m, the time step of the proposed method is about 3700 times as large as that of the FD method, and the total CPU time for the proposed method can be reduced to about 13% of the FD method.Table 2
The comparison of CPU time.
Δ
x
=
0.001 m
Δ
t
CPU time
FD method
1.33
×
1
0
-
7 s
237.9171 s
Proposed method
5
×
1
0
-
4 s
30.8813 sFigure 4
The numerical results at (a)p
1 and (b) p
2.
(a)
(b)


## 3.1. 1D Case

In order to validate the proposed method, a 1D structure constituted by two isotropic media separated by a thin layer (d
=
0.1 m, v
=
5.0 km/s) is analyzed. As shown in Figure 1, the thin layer is located at 60 m, and the total simulation length is 104 m. Displacement U
(
x
,
t
) at two points p
1 (38 m) and p
2 (78.1 m) is observed. The sound velocity in media A is 2.5 km/s, and the sound velocity in media B is 3.7 km/s. Nonuniform spatial grids are used. In the thin layer Δ
x
=
0.01 m; in other parts Δ
x
=
1 m.Figure 1
Configuration of a thin layer (d
=
0.1 m) embedded in two isotropic media.The external source is located atx
=
26 m and its expression is f
(
t
)
=
-
5.76
f
0
2
[
1
-
16
0.6
f
0
t
-
1
)
2
exp
⁡
[
-
8
(
0.6
f
0
t
-
1
)
2
], where f
0
=
100 Hz. Time support for the simulation is 0.16 s. To expand f
(
t
) in AH domain, we select N
=
52, σ
=
0.008, and T
f
=
0.5. Due to the limitation of the CFL condition, the time step is set as 2
×
1
0
-
6 s for the FD method under the condition of Δ
x
=
0.01 m. Time step is not involved in AH domain computation. However, for AH decomposition, we still need to specify a sample interval; the time step Δ
t
=
0.002 s is good enough. The absolute error R is defined as R
=
U
a
(
t
)
-
U
c
(
t
) for comparison of the two methods. U
c
(
t
) is the reference data obtained from the FD method with Δ
x
=
0.01 m, and U
a
(
t
) is the numerical result obtained from the FD method (Δ
x
=
0.1 m) and the proposed method (Δ
x
=
0.01 m). The numerical results at p
1 and their errors are shown in Figure 2(a). Figure 2(b) is the numerical results at p
2. From Figure 2, it can be seen that the proposed method has a good agreement with the FD method. The absolute error of the proposed method (Δ
x
=
0.01 m) is much lower than that of the FD method (Δ
x
=
0.1 m) at both p
1 and p
2.Figure 2
The numerical results and their errors at (a)p
1 and (b) p
2.
(a)
(b)Table1 shows the comparison of CPU time between the proposed method and the FD method. Under the condition of Δ
x
=
0.01 m, the time step of the proposed method is 1000 times as large as that of the FD method, and the total CPU time for the proposed method can be reduced to about 10% of the FD method.Table 1
The comparison of CPU time.
Δ
x = 0.01 m
Δ
t
CPU time
FD method
2
×
1
0
-
6 s
3.198 s
Proposed method
2
×
1
0
-
3 s
0.3127 s


## 3.2. 2D Case

In this section, we show results from a 2D numerical modeling by the FD method and the proposed method. The 2D model is established by two isotropic media separated by a thin layer (d
=
0.01 m, v
=
6.0 km/s), which is similar to the 1D case. As shown in Figure 3, the computational domain is 60 m × 90 m, and the thin layer is located at 72 m in x direction. Displacement U
(
x
,
y
,
t
) at two points p
1
(
36
m
,
30
m
) and p
2
(
81.01
m
,
30
m
) is observed. The sound velocity in media A is 3 km/s, and the sound velocity in media B is 4.7 km/s. Nonuniform spatial grids are used. In x direction, Δ
x
=
0.001 m in the thin layer, and Δ
x
=
3 m in other parts. In y direction, Δ
y
=
3 m. The external source is the same as the 1D case and is located at s
(
27
m
,
30
m
).Figure 3
Configuration of a thin layer (d
=
0.01 m) embedded in two isotropic media.Time support for the simulation is 0.06 s. The time step is set as1.33
×
1
0
-
7 s and 5
×
1
0
-
4 s for the FD method and the proposed method, respectively. Figure 4 shows the 2D modeling records at p
1 and p
2 for the multilayer isotropic media model. The waveform corresponding to the proposed method is in good agreement with that by the FD method, which demonstrates that the proposed method has good accuracy. Table 2 shows the comparison of CPU time between the proposed method and the FD method. Under the condition of Δ
x
=
0.001 m, the time step of the proposed method is about 3700 times as large as that of the FD method, and the total CPU time for the proposed method can be reduced to about 13% of the FD method.Table 2
The comparison of CPU time.
Δ
x
=
0.001 m
Δ
t
CPU time
FD method
1.33
×
1
0
-
7 s
237.9171 s
Proposed method
5
×
1
0
-
4 s
30.8813 sFigure 4
The numerical results at (a)p
1 and (b) p
2.
(a)
(b)


## 4. Conclusion

In this paper, an unconditionally stable method has been proposed for solving the 1D and 2D acoustic wave equation in time domain. By applying Associated Hermit orthogonal function expansion and Galerkin temporal testing procedure, the time variable can be eliminated from the calculations and can make the proposed method unconditionally stable. The numerical results show the proposed method is efficient and the accuracy is still guaranteed.



---

*Source: 381052-2015-08-11.xml*