# Generalized Finite Difference Time Domain Method and Its Application to Acoustics

Authors: Jianguo Wei; Song Wang; Qingzhi Hou; Jianwu Dang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A meshless generalized finite difference time domain (GFDTD) method is proposed and applied to transient acoustics to overcome difficulties due to use of grids or mesh. Inspired by the derivation of meshless particle methods, the generalized finite difference method (GFDM) is reformulated utilizing Taylor series expansion. It is in a way different from the conventional derivation of GFDM in which a weighted energy norm was minimized. The similarity and difference between GFDM and particle methods are hence conveniently examined. It is shown that GFDM has better performance than the modified smoothed particle method in approximating the first- and second-order derivatives of 1D and 2D functions. To solve acoustic wave propagation problems, GFDM is used to approximate the spatial derivatives and the leap-frog scheme is used for time integration.  By analog with FDTD, the whole algorithm is referred to as GFDTD. Examples in one- and two-dimensional domain with reflection and absorbing boundary conditions are solved and good agreements with the FDTD reference solutions are observed, even with irregular point distribution. The developed GFDTD method has advantages in solving wave propagation in domain with irregular and moving boundaries.


---

## Body

## 1. Introduction

Partial differential equations (PDEs) modeling problems in science and engineering, such as electromagnetics, acoustics, and hydrodynamics, are usually solved by numerical methods that discretize the computational domain with mesh or grids. Grid-based methods such as finite difference method (FDM), finite element method (FEM), and boundary element method (BEM) [1, 2] have had much achievements and still dominate the field of scientific computing. However, numerical difficulties originating from usage of grids often emerge. For complicated and irregular geometry, implementation of boundary conditions could be a big challenge for FDM. Generation of grids with high quality is not an easy task in FEM and BEM. Moreover, when free surface and moving boundary/interface have to be treated, the transformation of grids will turn the conventional grid-based methods into a difficult, time-consuming process. Numerical accuracy often degenerates and divergence problem occurs.In recent 20 years, to overcome numerical difficulties due to use of grids or mesh, meshless methods (MMs) based on different techniques have been proposed and widely used in many fields such as hydrodynamics [3], astrophysics [4], and solid mechanics [3, 5]. Among the MMs, generalized finite difference method (GFDM) is the one that evolved from traditional FDM [6, 7] and many different forms have been developed [8]. Benito and his coauthors made great contribution to its recent development [9–11]. For heat conduction problem, it has been compared with the element-free Garlerkin (EFG) method (one of the most used MMs in solid mechanics) and better performance has been observed [10]. Recently, GFDM was used to solve the wave equations [11] and Burgers’ equations [12] and simulate seismic wave propagation problems in heterogeneous media [13]. An application to the detonation shock dynamics [14] was also carried out. Nevertheless, few work on computational acoustics has been reported.For acoustic wave propagation problems, the concentration is on the ones in confined domain, for which grid-based methods like FDTD and TDFEM (time-domain finite-element methods) [15], are mostly used. However, moving boundary exists in many acoustic problems like sound wave propagation inside a deforming vocal tract. This problem is hardly solved by conventional grid-based methods and MMs provide a possibility. As one of the MMs, GFDM is extended to transient acoustics in this paper, which is helpful to solve wave propagation problems with moving boundary in the future.Inspired by the derivation of meshless particle methods, we firstly formulated the GFDM in a way different from the original one that minimizes an energy norm. Such that the relationship between GFDM and meshless particle methods like smoothed particle hydrodynamics (SPH) and its improvements can be conveniently examined. Comparison with the modified dmoothed particle hydrodynamics (MSPH) method, which has better performances than SPH and its corrections [16], shows higher approximation accuracy of the GFDM, especially at the boundary region. By analog with FDTD, a method referred to as generalized finite difference time domain (GFDTD) is proposed, in which GFDM is used to discretize the spatial operators and the leap-frog algorithm is used for time integration. To show its good performance and efficiency, the GFDTD method is applied to transient acoustics. Comparison with conventional FDTD solutions is presented and discussed.


## 2. Generalized Finite Difference Method (GFDM)

Other than conventional derivation of GFDM by minimizing an energy norm [10], a different derivation of GFDM is presented in this section. Taylor series expansion of f
(
x
,
y
) around point (
x
0
,
y
0
) remaining up to second-order terms yields(1)
f
≈
f
0
+
h
∂
f
0
∂
x
+
k
∂
f
0
∂
y
+
h
2
2
∂
2
f
0
∂
x
2
+
k
2
2
∂
2
f
0
∂
y
2
+
h
k
∂
2
f
0
∂
x
∂
y
,where f
=
f
(
x
,
y
), f
0
=
f
(
x
0
,
y
0
), h
=
x
-
x
0, and  k
=
y
-
y
0.By multiplying both sides of (1) with w
2
h and w
2
k (w is a weighting function with compact support) and integrating the resulted equations over the support domain Ω, we get two equations, and the following, as an example, is the result for w
2
h:(2)
∫
Ω
f
w
2
h
d
v
≈
∫
Ω
f
0
w
2
h
d
v
+
∫
Ω
∂
f
0
∂
x
w
2
h
2
d
v
+
∫
Ω
∂
f
0
∂
y
w
2
h
k
d
v
+
1
2
∫
Ω
∂
2
f
0
∂
x
2
w
2
h
3
d
v
+
1
2
∫
Ω
∂
2
f
0
∂
y
2
w
2
h
k
2
d
v
+
∫
Ω
∂
2
f
0
∂
x
∂
y
w
2
h
2
k
d
v
,where d
v is a volume measure.Repeating the same procedure withw
2
h
2
/
2, w
2
k
2
/
2, and w
2
h
k instead of w
2
h and w
2
k, we get other three equations and the following is the result for w
2
h
2
/
2: (3)
∫
Ω
f
w
2
h
2
2
d
v
≈
∫
Ω
f
0
w
2
h
2
2
d
v
+
∫
Ω
∂
f
0
∂
x
w
2
h
3
2
d
v
+
∫
Ω
∂
f
0
∂
y
w
2
h
2
k
2
d
v
+
1
2
∫
Ω
∂
2
f
0
∂
x
2
w
2
h
4
2
d
v
+
1
2
∫
Ω
∂
2
f
0
∂
y
2
w
2
h
2
k
2
2
d
v
+
∫
Ω
∂
2
f
0
∂
x
∂
y
w
2
h
3
k
2
d
v
.To approximate the integrations by Riemann sum, the volume of the support domain Ω is divided into N points with associated volumes d
v
i, (
i
=
1,2
,
…
,
N
). Equations (2), (3), and the other three constitute a system of five equations written in matrix form as(4)
A
P
D
f
P
=
b
P
,with(5)
A
P
=
∑
i
=
1
N
w
i
2
h
i
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
d
v
i
∑
i
=
1
N
w
i
2
h
i
3
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
2
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
d
v
i
∑
i
=
1
N
w
i
2
k
i
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
2
d
v
i
∑
i
=
1
N
w
i
2
k
i
3
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
3
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
2
d
v
i
∑
i
=
1
N
w
i
j
2
h
i
4
4
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
2
4
d
v
i
∑
i
=
1
N
w
i
2
h
i
3
h
i
2
d
v
i
∑
i
=
1
N
w
i
j
2
h
i
k
i
2
2
d
v
i
∑
i
=
1
N
w
i
2
k
i
3
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
2
4
d
v
i
∑
i
=
1
N
w
i
2
k
i
4
4
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
3
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
3
k
i
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
k
i
3
2
d
v
i
∑
i
=
1
N
w
i
2
h
i
2
k
i
2
d
v
i
,
D
f
P
=
∂
f
0
∂
x
∂
f
0
∂
y
∂
2
f
0
∂
x
2
∂
2
f
0
∂
y
2
∂
2
f
0
∂
x
∂
y
,
b
P
=
-
f
0
∑
i
=
1
N
w
i
2
h
i
d
v
i
+
∑
i
=
1
N
f
i
w
i
2
h
i
d
v
i
-
f
0
∑
i
=
1
N
w
i
2
k
i
d
v
i
+
∑
i
=
1
N
f
i
w
i
2
k
i
d
v
i
-
f
0
∑
i
=
1
N
w
i
2
h
i
2
2
d
v
i
+
∑
i
=
1
N
f
i
w
i
2
h
i
2
2
d
v
i
-
f
0
∑
i
=
1
N
w
i
2
k
i
2
2
d
v
i
+
∑
i
=
1
N
f
i
w
i
2
k
i
2
2
d
v
i
-
f
0
∑
i
=
1
N
w
i
2
h
i
k
i
d
v
i
+
∑
i
=
1
N
f
i
w
i
2
h
i
k
i
d
v
i
,where w
i
=
w
(
d
i
,
d
m
) with d
i
=
(
x
i
-
x
0
)
2
+
(
y
i
-
y
0
)
2, h
i
=
x
i
-
x
0, and k
i
=
y
i
-
y
0, and d
m is a measure of the support size.The conventional derivation of GFDM is presented in the appendix. It is clear that the difference between the conventional and the current derivation is not only the procedure but also the final form. The conventional derivation loses termd
v
i (see (A.5)). If all the points in the domain have the same volume, d
v
i at both sides of (4) will be cancelled, and the two final forms will be the same. However, d
v
i can hardly be the same when points are irregularly spaced. From this point of view, our derived final form is more general and takes point irregularity into account.


## 3. Modified Smoothed Particle Hydrodynamics (MSPH)

As a modification to SPH, the MSPH method improves the accuracy of the approximations especially at points near the boundary of the domain [16]. It uses Taylor series expansion of function f
(
x
,
y
) as in (1). Similar to the derivations of (2) and (3), but with different weight functions ∂
w
/
∂
x, ∂
w
/
∂
y, ∂
w
2
/
∂
x
2, ∂
w
2
/
∂
y
2, and ∂
2
w
/
∂
x
∂
y, the following equations, as examples, for ∂
w
/
∂
x and ∂
w
2
/
∂
x
2, are obtained:(6)
∫
Ω
f
∂
w
∂
x
d
v
=
∫
Ω
f
0
∂
w
∂
x
h
d
v
+
∂
f
0
∂
x
∫
Ω
∂
w
∂
x
h
d
v
+
∂
f
0
∂
y
∫
Ω
∂
w
∂
x
k
d
v
+
1
2
∂
2
f
0
∂
x
2
∫
Ω
∂
w
∂
x
h
2
d
v
+
1
2
∂
2
f
0
∂
y
2
∫
Ω
∂
w
∂
x
k
2
d
v
+
∂
2
f
0
∂
x
∂
y
∫
Ω
∂
w
∂
x
h
k
d
v
,
∫
Ω
f
∂
w
2
∂
x
2
d
v
=
∫
Ω
f
0
∂
w
2
∂
x
2
d
v
+
∂
f
∂
x
∫
Ω
∂
w
2
∂
x
2
h
d
v
+
∂
f
∂
y
∫
Ω
∂
w
2
∂
x
2
k
d
v
+
1
2
∂
2
f
0
∂
x
2
∫
Ω
∂
w
2
∂
x
2
h
2
d
v
+
1
2
∂
2
f
0
∂
y
2
∫
Ω
∂
w
2
∂
x
2
k
2
d
v
+
∂
2
f
0
∂
x
∂
y
∫
Ω
∂
w
2
∂
x
2
h
k
d
v
.Again the Riemann sum over the support domain Ω is used to approximate the integrations and a system of five equations is obtained as(7)
∑
i
=
1
N
∂
w
i
∂
x
h
i
d
v
i
∑
i
=
1
N
∂
w
i
∂
x
k
i
d
v
i
∑
i
=
1
N
∂
w
i
∂
x
h
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
∂
x
k
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
∂
x
h
i
k
i
d
v
i
∑
i
=
1
N
∂
w
i
∂
y
h
i
d
v
i
∑
i
=
1
N
∂
w
i
∂
y
k
i
d
v
i
∑
i
=
1
N
∂
w
i
∂
y
h
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
∂
y
k
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
∂
y
h
i
k
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
2
h
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
2
k
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
2
h
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
2
k
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
2
h
i
k
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
y
2
h
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
y
2
k
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
y
2
h
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
y
2
k
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
y
2
h
i
k
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
∂
y
h
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
∂
y
k
i
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
∂
y
h
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
∂
y
k
i
2
2
d
v
i
∑
i
=
1
N
∂
w
i
2
∂
x
∂
y
h
i
k
i
d
v
i
∂
f
0
∂
x
∂
f
0
∂
y
∂
2
f
0
∂
x
2
∂
2
f
0
∂
y
2
∂
2
f
0
∂
x
∂
y
=
-
f
0
∑
i
=
1
N
∂
w
i
∂
x
d
v
i
+
∑
i
=
1
N
f
i
∂
w
i
∂
x
d
v
i
-
f
0
∑
i
=
1
N
∂
w
i
∂
y
d
v
i
+
∑
i
=
1
N
f
i
∂
w
i
∂
y
d
v
i
-
f
0
∑
i
=
1
N
∂
w
i
2
∂
x
2
d
v
i
+
∑
i
=
1
N
f
i
∂
w
i
2
∂
x
2
d
v
i
-
f
0
∑
i
=
1
N
∂
w
i
2
∂
y
2
d
v
i
+
∑
i
=
1
N
f
i
∂
w
i
2
∂
y
2
d
v
i
-
f
0
∑
i
=
1
N
∂
w
i
2
∂
x
∂
y
d
v
i
+
∑
i
=
1
N
f
i
∂
w
i
2
∂
x
∂
y
d
v
i
.Compared with formula (4), the only difference is the terms multiplied to both sides of (1). In GFDM, w
2
h, w
2
k, w
2
h
2
/
2, w
2
k
2
/
2, and w
2
h
k are used instead of ∂
w
/
∂
x, ∂
w
/
∂
y, ∂
w
2
/
∂
x
2, ∂
w
2
/
∂
y
2, and ∂
2
w
/
∂
x
∂
y in MSPH. As a result, GFDM avoids computing the derivatives of the weight function and hence saves computational efforts and leads to more choice of the weight function.


## 4. Numerical Tests for Approximation of Derivatives

In previous sections the deviation of GFDM and MSPH is presented. In this section, to compare the performance of the two methods, they are used to approximate the derivatives of certain 1D and 2D functions. For the convenience of evaluation, a global error measure is defined as follows:(8)
Error
u
=
1
u
max
⁡
1
N
∑
i
=
1
N
u
i
e
-
u
i
n
2
,where u can be ∂
f
/
∂
x, ∂
f
/
∂
y, ∂
2
f
/
∂
x
2, and ∂
2
f
/
∂
y
2 and the superscripts (
e
) and (
n
) refer to the exact and numerical solutions, respectively.The quartic spline function is used as the weight functionw
i: (9)
w
i
(
d
)
=
1
-
6
d
d
m
2
+
8
d
d
m
3
-
3
d
d
m
4
,
d
≤
d
m
,
0
,
d
>
d
m
,where d
m is the kernel radius taken as 2.1
Δ
x (Δ
x is the space interval) which is usually used in meshless methods.


### 4.1. One-Dimensional Case

Consider the following function:(10)
f
(
x
)
=
x
-
0.5
4
,
x
∈
0,1
.Figure 1 shows the first- and second-order derivatives estimated by GFDM and MSPH and the exact results when the domain is discretized into 21 equally spaced points. It is seen that GFDM has better performance in both derivatives especially for the points near boundaries. When the number of points increases to 51, the results are similar as exhibited in Figure 2. Error analysis shown in Table 1 indicates that GFDM has higher accuracy. With increasing number of points, the global error decreases.Table 1
Approximation errors in the derivatives of the functionf
(
x
)
=
(
x
-
0.5
)
4.
Error (%)
∂
f
/
∂
x
∂
2
f
/
∂
x
2
GFDM with 21 points
0.77
5.81
MSPH with 21 points
0.77
6.02
GFDM with 51 points
0.11
1.55
MSPH with 51 points
0.11
1.61Figure 1
Estimates of the first derivative (a) and second derivative (b) of functionf
(
x
)
=
(
x
-
0.5
)
4 with 21 equally spaced points on [
0
,
1
].
(a)
(b)Figure 2
Estimates of the first derivative (a) and second derivative (b) of functionf
(
x
)
=
(
x
-
0.5
)
4 with 51 equally spaced points on [
0
,
1
].
(a)
(b)


### 4.2. Two-Dimensional Case

For the function(11)
f
(
x
,
y
)
=
sin
π
x
sin
π
y
,
x
,
y
∈
[
0,1
]
×
0,1
,its first- and second-order derivatives together with estimations by GFDM and MSPH are shown in Figure 3. In each direction 21 points are employed. As expected, GFDM has higher approximation accuracy than MSPH for both first- and second-order derivatives as shown in Table 2.Table 2
Approximation errors in the derivatives of the functionf
(
x
,
y
)
=
sin
⁡
⁡
π
x
sin
⁡
⁡
π
y.
Error (%)
∂
f
/
∂
x
∂
f
/
∂
y
∂
2
f
/
∂
x
2
∂
2
f
/
∂
y
2
GFDM
0.29
0.29
2.87
2.87
MSPH
0.47
0.47
4.71
4.71Figure 3
Exact (top) and estimated (a)∂
f
/
∂
x and (b) ∂
2
f
/
∂
x
2 of function f
(
x
,
y
)
=
sin
π
x
sin
π
y by GFDM (middle) and MSPH (bottom).
(a)
(b)


## 4.1. One-Dimensional Case

Consider the following function:(10)
f
(
x
)
=
x
-
0.5
4
,
x
∈
0,1
.Figure 1 shows the first- and second-order derivatives estimated by GFDM and MSPH and the exact results when the domain is discretized into 21 equally spaced points. It is seen that GFDM has better performance in both derivatives especially for the points near boundaries. When the number of points increases to 51, the results are similar as exhibited in Figure 2. Error analysis shown in Table 1 indicates that GFDM has higher accuracy. With increasing number of points, the global error decreases.Table 1
Approximation errors in the derivatives of the functionf
(
x
)
=
(
x
-
0.5
)
4.
Error (%)
∂
f
/
∂
x
∂
2
f
/
∂
x
2
GFDM with 21 points
0.77
5.81
MSPH with 21 points
0.77
6.02
GFDM with 51 points
0.11
1.55
MSPH with 51 points
0.11
1.61Figure 1
Estimates of the first derivative (a) and second derivative (b) of functionf
(
x
)
=
(
x
-
0.5
)
4 with 21 equally spaced points on [
0
,
1
].
(a)
(b)Figure 2
Estimates of the first derivative (a) and second derivative (b) of functionf
(
x
)
=
(
x
-
0.5
)
4 with 51 equally spaced points on [
0
,
1
].
(a)
(b)


## 4.2. Two-Dimensional Case

For the function(11)
f
(
x
,
y
)
=
sin
π
x
sin
π
y
,
x
,
y
∈
[
0,1
]
×
0,1
,its first- and second-order derivatives together with estimations by GFDM and MSPH are shown in Figure 3. In each direction 21 points are employed. As expected, GFDM has higher approximation accuracy than MSPH for both first- and second-order derivatives as shown in Table 2.Table 2
Approximation errors in the derivatives of the functionf
(
x
,
y
)
=
sin
⁡
⁡
π
x
sin
⁡
⁡
π
y.
Error (%)
∂
f
/
∂
x
∂
f
/
∂
y
∂
2
f
/
∂
x
2
∂
2
f
/
∂
y
2
GFDM
0.29
0.29
2.87
2.87
MSPH
0.47
0.47
4.71
4.71Figure 3
Exact (top) and estimated (a)∂
f
/
∂
x and (b) ∂
2
f
/
∂
x
2 of function f
(
x
,
y
)
=
sin
π
x
sin
π
y by GFDM (middle) and MSPH (bottom).
(a)
(b)


## 5. Generalized Finite Difference Time Domain Method for Computational Acoustics

For computational acoustics, the mostly used approach is the FDTD method, which was originally designed for the simulation of electromagnetics [1, 2]. As a finite difference scheme, its applicability to complex problems suffers from aforementioned difficulties, for which the generalized finite difference can be a good alternative. In this section, together with the basics of computational acoustics, a meshless method is proposed, in which GFDM is used to discretize the spatial derivatives and the leap-frog algorithm is used to discretize the temporal derivatives. By analog with FDTD, it is referred to as generalized finite difference time domain (GFDTD) method and is expected to have advantages due to its meshless property.The governing equations for acoustic wave propagation problems are(12)
ρ
0
∂
v
∂
t
=
-
∇
p
,
1
c
0
2
∂
p
∂
t
=
-
ρ
0
∇
·
v
,where p is pressure, v is particle velocity, ρ
0 is the density of the medium, and c
0 is the speed of sound.


### 5.1. Spatial Derivative Approximations by GFDM

The spatial derivatives on the right-hand side of (12) are approximated by GFDM. By solving (4) we get the approximations of ∂
f
/
∂
x and ∂
f
/
∂
y. That is, the derivatives of variable f at point x
0
,
y
0 can be approximated by function values at points inside the support domain centered at x
0
,
y
0 as(13)
∂
f
0
∂
x
=
-
m
0
f
0
+
∑
i
=
1
N
f
i
m
i
,
∂
f
0
∂
y
=
-
η
0
f
0
+
∑
i
=
1
N
f
i
η
i
,where m
0 and η
0 are the difference coefficients for center point and m
i and η
i are coefficients for other points in the support domain. As GFDM can reproduce constant functions [4], we have(14)
m
0
=
∑
i
=
1
N
m
i
,
η
0
=
∑
i
=
1
N
η
i
.By taking both p and v in (12) as f, the approximated spatial operators in (12) are accordingly obtained.


### 5.2. Explicit Leap-Frog Scheme in GFDTD

Generally, the temporal derivatives on the left-hand side of (12) can be integrated by any time marching algorithms. Inspired by the conventional FDTD method, the second-order accurate explicit leap-frog scheme is used herein, in which two variables p and v are alternatively calculated. The velocity is computed at the half time step and the pressure is calculated at the integer time step [2]. After temporal approximations the semi-discretization of (12) becomes(15)
ρ
0
v
n
+
1
/
2
-
v
n
-
1
/
2
Δ
t
=
-
∇
p
n
,
1
c
0
2
p
n
+
1
-
p
n
Δ
t
=
-
ρ
0
∇
·
v
n
+
1
/
2
,where superscript n represents the time step. The leap-frog scheme is conditionally stable and the time step Δ
t should satisfy the Courant-Friedrichs-Lewy (CFL) condition; that is, Δ
t
≤
Δ
x
/
(
dim
⁡
·
c
0
), where dim
⁡ is the dimension of the problem.Substituting (13) into the right-hand side of (15) as the approximation of the first-order spatial derivatives, the full discretized system of equations becomes(16)
v
0
n
+
1
/
2
=
v
0
n
-
1
/
2
+
Δ
t
ρ
0
·
m
0
p
0
n
-
∑
i
=
1
N
m
i
p
i
n
,
η
0
p
0
n
-
∑
i
=
1
N
η
i
p
i
n
T
,
p
0
n
+
1
=
p
0
n
+
Δ
t
ρ
0
c
0
2
·
m
0
u
0
n
+
1
/
2
-
∑
i
=
1
N
m
i
u
i
n
+
1
/
2
+
η
0
v
0
n
+
1
/
2
-
∑
i
=
1
N
η
i
v
i
n
+
1
/
2
,where the particle velocity is a 2D vector v
=
u
,
v
T.By analog with FDTD, the full discretization scheme is referred to as generalized finite difference time domain (GFDTD) method.


## 5.1. Spatial Derivative Approximations by GFDM

The spatial derivatives on the right-hand side of (12) are approximated by GFDM. By solving (4) we get the approximations of ∂
f
/
∂
x and ∂
f
/
∂
y. That is, the derivatives of variable f at point x
0
,
y
0 can be approximated by function values at points inside the support domain centered at x
0
,
y
0 as(13)
∂
f
0
∂
x
=
-
m
0
f
0
+
∑
i
=
1
N
f
i
m
i
,
∂
f
0
∂
y
=
-
η
0
f
0
+
∑
i
=
1
N
f
i
η
i
,where m
0 and η
0 are the difference coefficients for center point and m
i and η
i are coefficients for other points in the support domain. As GFDM can reproduce constant functions [4], we have(14)
m
0
=
∑
i
=
1
N
m
i
,
η
0
=
∑
i
=
1
N
η
i
.By taking both p and v in (12) as f, the approximated spatial operators in (12) are accordingly obtained.


## 5.2. Explicit Leap-Frog Scheme in GFDTD

Generally, the temporal derivatives on the left-hand side of (12) can be integrated by any time marching algorithms. Inspired by the conventional FDTD method, the second-order accurate explicit leap-frog scheme is used herein, in which two variables p and v are alternatively calculated. The velocity is computed at the half time step and the pressure is calculated at the integer time step [2]. After temporal approximations the semi-discretization of (12) becomes(15)
ρ
0
v
n
+
1
/
2
-
v
n
-
1
/
2
Δ
t
=
-
∇
p
n
,
1
c
0
2
p
n
+
1
-
p
n
Δ
t
=
-
ρ
0
∇
·
v
n
+
1
/
2
,where superscript n represents the time step. The leap-frog scheme is conditionally stable and the time step Δ
t should satisfy the Courant-Friedrichs-Lewy (CFL) condition; that is, Δ
t
≤
Δ
x
/
(
dim
⁡
·
c
0
), where dim
⁡ is the dimension of the problem.Substituting (13) into the right-hand side of (15) as the approximation of the first-order spatial derivatives, the full discretized system of equations becomes(16)
v
0
n
+
1
/
2
=
v
0
n
-
1
/
2
+
Δ
t
ρ
0
·
m
0
p
0
n
-
∑
i
=
1
N
m
i
p
i
n
,
η
0
p
0
n
-
∑
i
=
1
N
η
i
p
i
n
T
,
p
0
n
+
1
=
p
0
n
+
Δ
t
ρ
0
c
0
2
·
m
0
u
0
n
+
1
/
2
-
∑
i
=
1
N
m
i
u
i
n
+
1
/
2
+
η
0
v
0
n
+
1
/
2
-
∑
i
=
1
N
η
i
v
i
n
+
1
/
2
,where the particle velocity is a 2D vector v
=
u
,
v
T.By analog with FDTD, the full discretization scheme is referred to as generalized finite difference time domain (GFDTD) method.


## 6. Numerical Results

To validate the proposed GFDTD method, it is applied to one- and two-dimensional wave propagation problems. Three cases are presented. The first two examine the acoustic wave propagation in one- and two-dimensional domain, respectively, and the third one is a real case with different types of boundary conditions. All the cases use (9) as the weight function. The other parameters are ρ
0
=
1
kg
/
m
3 and c
0
=
346.4 m/s and the time interval Δ
t is set to be 1 μs to satisfy the CFL condition. FDTD solutions are chosen as the reference.


### 6.1. One-Dimensional Case

In this case 501 points are equally spaced in the domain[
0
,
1
]. In the middle of it, there is a wave source in Gaussian pulse form:(17)
gp
(
x
)
=
e
-
25
x
-
0.5
.As shown in Figure4, the simulated results at two time levels t
=
250 μs and 500 μs have good agreement with the FDTD solutions. Compared with FDTD, the relative errors concerning the pressure P
FDTD and P
GFDTD are(18)
Error
250
=
P
FDTD
-
P
GFDTD
2
P
FDTD
2
=
0.267
·
1
0
-
2
,
Error
500
=
P
FDTD
-
P
GFDTD
2
P
FDTD
2
=
0.519
·
1
0
-
2
.Figure 4
Gaussian wave propagation at two time points: (a)t
=
250 μs and (b) 500 μs.
(a)
(b)


### 6.2. Two-Dimensional Case

The Gaussian wave propagation in two-dimensional domain is simulated in this section. The length of the square domain is 0.1 m and 101 points are uniformly distributed in each direction. The wave source starts from the middle of the domain. Figure5 compares the solutions of FDTD and GFDTD at t
=
20 μs. The results along y
=
0.05 m are shown in Figure 5(c). Again, good agreement is observed and the relative error is less than 2%.Figure 5
Gaussian wave propagation in a square domain at timet
=
20 μs. (a) The FDTD result, (b) the GFDTD result, and (c) the comparison along y
=
0.05 m.
(a)
(b)
(c)To show the advantage of the proposed GFDTD over the conventional FDTD, irregular point distribution is examined. All the points used above are allowed to have±
10
% perturbation around their original locations to make the distribution irregular, part of which is shown in Figure 6(a) and the result is shown in Figure 6(b). The comparison of the result along y
=
0.05 m shown in Figure 6(c) indicates that, with irregular distribution of computational points, the Gaussian wave propagates as well as before.Figure 6
Gaussian wave propagation in a square domain with irregular point distribution. (a) Part of the irregular distribution, (b) simulated results, and (c) comparison with regular point distribution after20 μs.
(a)
(b)
(c)In the GFDTD simulation of wave propagation with irregular point distribution, the volume associated to each point had better to be considered as analyzed at the end of Section2. In 2D case, the volume associated with a given point is the area that the point dominates. Here we use Delaunay triangulation and Voronoi diagram [17] to calculate the area and the results are shown in Figure 7. The volume of each point is shown in Figure 7(a). Due to the designed ±
10
% perturbation, the volume fluctuates around 1
-
6 m2. Based on Delaunay triangulation and Voronoi diagram, the area associated with a red point is indicated by the area surrounded by the blue lines as shown in Figure 7(b). To compare the results with regular and irregular point distributions, the values on regularly distributed points have to be firstly interpolated by the calculated results with irregularly spaced points. The third-order accurate cubic interpolation method [17] is used herein. The relative error is (19)
Error
=
P
iriegular
-
P
regular
2
P
regular
2
=
1.58
·
1
0
-
2
.The relatively small error is due to the gentle irregular point distribution. That is, if the point irregularity is small, the effect of d
v
i is negligible. This is consistent with Benito’s results [9–11].Figure 7
(a) Area associated with irregular distributed points and (b) Voronoi diagram.
(a)
(b)


### 6.3. Two-Dimensional Case with Effect of Boundary Conditions

In Sections6.1 and 6.2, wave propagation inside a domain is simulated. In this section, to show the effect of boundaries, which is of high importance in transient acoustics, sound wave propagation in a rectangular tube is studied. Two different kinds of boundary conditions including reflection and absorbing boundary are considered. At the left edge of the computational domain there is a Gaussian pulse as the source term. The upper and bottom boundaries are reflection layers and the second-order Mur’s absorbing boundary condition [18] exists at the right side boundary. In this case, the same ρ
0 and c
0 are used as before and the time interval Δ
t is 1 μs. Inside the computational domain 64
×
100 points with spatial interval Δ
x
=
Δ
y
=
1 mm are evenly spaced.


#### 6.3.1. Source Term

The left is a wave source with pressure given by the Gaussian pulse:(20)
gp
(
t
)
=
e
-
t
-
T
/
0.29
T
2
,where T
=
0.646
/
f
0 and f
0
=
10 KHz.


#### 6.3.2. Reflection Boundary Condition

To simulate reflections at the upper and bottom wall boundaries, the model proposed by Yokota et al. [19, 20] and widely used in room acoustics is employed herein. In this model, the normal component of particle velocity and the pressure of the points on the boundary are supposed to satisfy the following condition: (21)
v
norm
=
p
Z
norm
,where Z
norm is the normal acoustic impedance on the boundary given by(22)
Z
norm
=
ρ
0
c
0
1
+
1
-
α
norm
1
-
1
-
α
norm
.Here the normal sound absorption coefficient α
norm is taken as 0.2 as in [19].


#### 6.3.3. Absorbing Boundary Condition

At the right boundary, second-order Mur’s absorbing boundary condition [18] is applied:(23)
1
c
0
∂
2
p
∂
x
∂
t
+
1
c
0
2
∂
2
p
∂
t
2
+
1
2
∂
2
p
∂
y
2
=
0
.By applying (12) to (23) and performing time integration, (23) degenerates to(24)
∂
p
∂
x
+
1
c
0
∂
p
∂
t
-
c
0
ρ
0
2
∂
u
∂
y
=
0
.When GFDTD is used, the discrete form of (24) is obtained as (25)
p
0
n
+
1
-
p
0
n
Δ
t
=
c
0
2
ρ
0
2
η
0
u
0
n
+
1
/
2
-
∑
i
=
1
N
η
i
u
i
n
+
1
/
2
-
c
0
m
0
p
0
n
-
∑
i
=
1
N
m
i
p
i
n
.


#### 6.3.4. Results

After applying the source term and the two boundary conditions into our case, the wave is considered to be propagating from left to right inside a tube and gets absorbed at the end of it. Figure8 shows the simulated results after 200 μs with a color map image that clearly depicts the pressure distribution. In Figure 9, the results after 350 μs are depicted and the absorbing boundary at the right edge leads to no reflection.Figure 8
Sound wave propagation in a two-dimension tube. (a) FDTD results, (b) GFDTD results, and (c) the comparison alongy
=
0.034m at t
=
200 μs.
(a)
(b)
(c)Figure 9
Sound wave propagation in a two-dimension tube. (a) FDTD results, (b) GFDTD results, and (c) the comparison alongy
=
0.034 m at t
=
350 μs.
(a)
(b)
(c)


## 6.1. One-Dimensional Case

In this case 501 points are equally spaced in the domain[
0
,
1
]. In the middle of it, there is a wave source in Gaussian pulse form:(17)
gp
(
x
)
=
e
-
25
x
-
0.5
.As shown in Figure4, the simulated results at two time levels t
=
250 μs and 500 μs have good agreement with the FDTD solutions. Compared with FDTD, the relative errors concerning the pressure P
FDTD and P
GFDTD are(18)
Error
250
=
P
FDTD
-
P
GFDTD
2
P
FDTD
2
=
0.267
·
1
0
-
2
,
Error
500
=
P
FDTD
-
P
GFDTD
2
P
FDTD
2
=
0.519
·
1
0
-
2
.Figure 4
Gaussian wave propagation at two time points: (a)t
=
250 μs and (b) 500 μs.
(a)
(b)


## 6.2. Two-Dimensional Case

The Gaussian wave propagation in two-dimensional domain is simulated in this section. The length of the square domain is 0.1 m and 101 points are uniformly distributed in each direction. The wave source starts from the middle of the domain. Figure5 compares the solutions of FDTD and GFDTD at t
=
20 μs. The results along y
=
0.05 m are shown in Figure 5(c). Again, good agreement is observed and the relative error is less than 2%.Figure 5
Gaussian wave propagation in a square domain at timet
=
20 μs. (a) The FDTD result, (b) the GFDTD result, and (c) the comparison along y
=
0.05 m.
(a)
(b)
(c)To show the advantage of the proposed GFDTD over the conventional FDTD, irregular point distribution is examined. All the points used above are allowed to have±
10
% perturbation around their original locations to make the distribution irregular, part of which is shown in Figure 6(a) and the result is shown in Figure 6(b). The comparison of the result along y
=
0.05 m shown in Figure 6(c) indicates that, with irregular distribution of computational points, the Gaussian wave propagates as well as before.Figure 6
Gaussian wave propagation in a square domain with irregular point distribution. (a) Part of the irregular distribution, (b) simulated results, and (c) comparison with regular point distribution after20 μs.
(a)
(b)
(c)In the GFDTD simulation of wave propagation with irregular point distribution, the volume associated to each point had better to be considered as analyzed at the end of Section2. In 2D case, the volume associated with a given point is the area that the point dominates. Here we use Delaunay triangulation and Voronoi diagram [17] to calculate the area and the results are shown in Figure 7. The volume of each point is shown in Figure 7(a). Due to the designed ±
10
% perturbation, the volume fluctuates around 1
-
6 m2. Based on Delaunay triangulation and Voronoi diagram, the area associated with a red point is indicated by the area surrounded by the blue lines as shown in Figure 7(b). To compare the results with regular and irregular point distributions, the values on regularly distributed points have to be firstly interpolated by the calculated results with irregularly spaced points. The third-order accurate cubic interpolation method [17] is used herein. The relative error is (19)
Error
=
P
iriegular
-
P
regular
2
P
regular
2
=
1.58
·
1
0
-
2
.The relatively small error is due to the gentle irregular point distribution. That is, if the point irregularity is small, the effect of d
v
i is negligible. This is consistent with Benito’s results [9–11].Figure 7
(a) Area associated with irregular distributed points and (b) Voronoi diagram.
(a)
(b)


## 6.3. Two-Dimensional Case with Effect of Boundary Conditions

In Sections6.1 and 6.2, wave propagation inside a domain is simulated. In this section, to show the effect of boundaries, which is of high importance in transient acoustics, sound wave propagation in a rectangular tube is studied. Two different kinds of boundary conditions including reflection and absorbing boundary are considered. At the left edge of the computational domain there is a Gaussian pulse as the source term. The upper and bottom boundaries are reflection layers and the second-order Mur’s absorbing boundary condition [18] exists at the right side boundary. In this case, the same ρ
0 and c
0 are used as before and the time interval Δ
t is 1 μs. Inside the computational domain 64
×
100 points with spatial interval Δ
x
=
Δ
y
=
1 mm are evenly spaced.


### 6.3.1. Source Term

The left is a wave source with pressure given by the Gaussian pulse:(20)
gp
(
t
)
=
e
-
t
-
T
/
0.29
T
2
,where T
=
0.646
/
f
0 and f
0
=
10 KHz.


### 6.3.2. Reflection Boundary Condition

To simulate reflections at the upper and bottom wall boundaries, the model proposed by Yokota et al. [19, 20] and widely used in room acoustics is employed herein. In this model, the normal component of particle velocity and the pressure of the points on the boundary are supposed to satisfy the following condition: (21)
v
norm
=
p
Z
norm
,where Z
norm is the normal acoustic impedance on the boundary given by(22)
Z
norm
=
ρ
0
c
0
1
+
1
-
α
norm
1
-
1
-
α
norm
.Here the normal sound absorption coefficient α
norm is taken as 0.2 as in [19].


### 6.3.3. Absorbing Boundary Condition

At the right boundary, second-order Mur’s absorbing boundary condition [18] is applied:(23)
1
c
0
∂
2
p
∂
x
∂
t
+
1
c
0
2
∂
2
p
∂
t
2
+
1
2
∂
2
p
∂
y
2
=
0
.By applying (12) to (23) and performing time integration, (23) degenerates to(24)
∂
p
∂
x
+
1
c
0
∂
p
∂
t
-
c
0
ρ
0
2
∂
u
∂
y
=
0
.When GFDTD is used, the discrete form of (24) is obtained as (25)
p
0
n
+
1
-
p
0
n
Δ
t
=
c
0
2
ρ
0
2
η
0
u
0
n
+
1
/
2
-
∑
i
=
1
N
η
i
u
i
n
+
1
/
2
-
c
0
m
0
p
0
n
-
∑
i
=
1
N
m
i
p
i
n
.


### 6.3.4. Results

After applying the source term and the two boundary conditions into our case, the wave is considered to be propagating from left to right inside a tube and gets absorbed at the end of it. Figure8 shows the simulated results after 200 μs with a color map image that clearly depicts the pressure distribution. In Figure 9, the results after 350 μs are depicted and the absorbing boundary at the right edge leads to no reflection.Figure 8
Sound wave propagation in a two-dimension tube. (a) FDTD results, (b) GFDTD results, and (c) the comparison alongy
=
0.034m at t
=
200 μs.
(a)
(b)
(c)Figure 9
Sound wave propagation in a two-dimension tube. (a) FDTD results, (b) GFDTD results, and (c) the comparison alongy
=
0.034 m at t
=
350 μs.
(a)
(b)
(c)


## 6.3.1. Source Term

The left is a wave source with pressure given by the Gaussian pulse:(20)
gp
(
t
)
=
e
-
t
-
T
/
0.29
T
2
,where T
=
0.646
/
f
0 and f
0
=
10 KHz.


## 6.3.2. Reflection Boundary Condition

To simulate reflections at the upper and bottom wall boundaries, the model proposed by Yokota et al. [19, 20] and widely used in room acoustics is employed herein. In this model, the normal component of particle velocity and the pressure of the points on the boundary are supposed to satisfy the following condition: (21)
v
norm
=
p
Z
norm
,where Z
norm is the normal acoustic impedance on the boundary given by(22)
Z
norm
=
ρ
0
c
0
1
+
1
-
α
norm
1
-
1
-
α
norm
.Here the normal sound absorption coefficient α
norm is taken as 0.2 as in [19].


## 6.3.3. Absorbing Boundary Condition

At the right boundary, second-order Mur’s absorbing boundary condition [18] is applied:(23)
1
c
0
∂
2
p
∂
x
∂
t
+
1
c
0
2
∂
2
p
∂
t
2
+
1
2
∂
2
p
∂
y
2
=
0
.By applying (12) to (23) and performing time integration, (23) degenerates to(24)
∂
p
∂
x
+
1
c
0
∂
p
∂
t
-
c
0
ρ
0
2
∂
u
∂
y
=
0
.When GFDTD is used, the discrete form of (24) is obtained as (25)
p
0
n
+
1
-
p
0
n
Δ
t
=
c
0
2
ρ
0
2
η
0
u
0
n
+
1
/
2
-
∑
i
=
1
N
η
i
u
i
n
+
1
/
2
-
c
0
m
0
p
0
n
-
∑
i
=
1
N
m
i
p
i
n
.


## 6.3.4. Results

After applying the source term and the two boundary conditions into our case, the wave is considered to be propagating from left to right inside a tube and gets absorbed at the end of it. Figure8 shows the simulated results after 200 μs with a color map image that clearly depicts the pressure distribution. In Figure 9, the results after 350 μs are depicted and the absorbing boundary at the right edge leads to no reflection.Figure 8
Sound wave propagation in a two-dimension tube. (a) FDTD results, (b) GFDTD results, and (c) the comparison alongy
=
0.034m at t
=
200 μs.
(a)
(b)
(c)Figure 9
Sound wave propagation in a two-dimension tube. (a) FDTD results, (b) GFDTD results, and (c) the comparison alongy
=
0.034 m at t
=
350 μs.
(a)
(b)
(c)


## 7. Conclusion

A new derivation of the generalized finite difference method (GFDM) with Taylor series expansion generates the same formulation as its conventional derivation and clearly demonstrates its relationship with meshless particle methods. GFDM has better performance in derivative approximations than the particle methods. The proposed generalized finite difference time domain (GFDTD) method has been successfully applied to one- and two-dimensional acoustic wave propagation problems with reflection and absorbing boundary conditions. The numerical results are in line with the FDTD reference solutions even with irregular point distribution. The GFDTD method has high potentials in solving transient acoustic problems with moving boundaries, which deserves further studies.



---

*Source: 640305-2015-06-23.xml*