Generalized Finite Difference Time Domain Method and Its Application to Acoustics

Ameshless 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 firstand 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 oneand 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.


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][10][11].For heat conduction problem, it has been compared with the elementfree 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 gridbased methods like FDTD and TDFEM (time-domain finiteelement 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 leapfrog 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.

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 (, ) around point ( 0 ,  0 ) remaining up to second-order terms yields where  = (, ),  0 = ( 0 ,  0 ), ℎ =  −  0 , and  =  −  0 .By multiplying both sides of (1) with  2 ℎ and  2  ( 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  2 ℎ: where dV is a volume measure.
The conventional derivation of GFDM is presented in 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 term dV  (see (A.5)).If all the points in the domain have the same volume, dV  at both sides of (4) will be cancelled, and the two final forms will be the same.However, dV  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.

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 (, ) as in (1).Similar to the derivations of ( 2) and ( 3), but with different weight functions /, /, Again the Riemann sum over the support domain Ω is used to approximate the integrations and a system of five equations is obtained as ) .
Compared with formula (4), the only difference is the terms multiplied to both sides of (1).In GFDM,  2 ℎ,

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: where  can be /, /,  2 / 2 , and  2 / 2 and the superscripts () and () refer to the exact and numerical solutions, respectively.The quartic spline function is used as the weight function   : where   is the kernel radius taken as 2.1Δ (Δ is the space interval) which is usually used in meshless methods.
4.1.One-Dimensional Case.Consider the following function: 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.

Two-Dimensional Case.
For the function 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 firstand second-order derivatives as shown in Table 2.

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  where  is pressure, k is particle velocity,  0 is the density of the medium, and  0 is the speed of sound.

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 / and /.That is, the derivatives of variable  at point ( 0 ,  0 ) can be approximated by function values at points inside the support domain centered at ( 0 ,  0 ) as where  0 and  0 are the difference coefficients for center point and   and   are coefficients for other points in the support domain.As GFDM can reproduce constant functions [4], we have By taking both  and k in (12) as , the approximated spatial operators in ( 12) are accordingly obtained.

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  and k 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 where superscript  represents the time step.The leap-frog scheme is conditionally stable and the time step Δ should satisfy the Courant-Friedrichs-Lewy (CFL) condition; that is, Δ ≤ Δ/( √ dim ⋅  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 where the particle velocity is a 2D vector k = [, V]  .By analog with FDTD, the full discretization scheme is referred to as generalized finite difference time domain (GFDTD) method.

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  0 = 346.4m/s and the time interval Δ 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: As shown in Figure 4, the simulated results at two time levels  = 250 s and 500 s have good agreement with the FDTD solutions.Compared with FDTD, the relative errors concerning the pressure  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  the result is shown in Figure 6(b).The comparison of the result along  = 0.05 m shown in Figure 6(c) indicates that, with irregular distribution of computational points, the Gaussian wave propagates as well as before.
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 Section 2. 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 perturbation, the volume fluctuates around 1 −6 m 2 .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] The relatively small error is due to the gentle irregular point distribution.That is, if the point irregularity is small, the effect of dV  is negligible.This is consistent with Benito's results [9][10][11].
6.3.Two-Dimensional Case with Effect of Boundary Conditions.In Sections 6.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  0 are used as before and the time interval Δ is 1 s.Inside the computational domain 64 × 100 points with spatial interval Δ = Δ = 1 mm are evenly spaced.

Source Term.
The left is a wave source with pressure given by the Gaussian pulse: where  = 0.646/ 0 and  0 = 10 KHz.

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: where  norm is the normal acoustic impedance on the boundary given by Here the normal sound absorption coefficient  norm is taken as 0.2 as in [19].

Absorbing Boundary Condition.
At the right boundary, second-order Mur's absorbing boundary condition [18] is applied: By applying ( 12) to (23) and performing time integration, (23) degenerates to When GFDTD is used, the discrete form of ( 24) is obtained as (25) 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.Figure 8 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.

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 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.

2 , 2 . ( 18 ) 6 . 2 .
FDTD and  GFDTD are Error 250 =       FDTD −  GFDTD     2      FDTD    2 = 0.267 ⋅ 10 −Error 500 =       FDTD −  GFDTD     2      FDTD    2 = 0.519 ⋅ 10 −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.Figure 5 compares the solutions of FDTD and GFDTD at  = 20 s.The results along  = 0.05 m are shown in Figure 5(c).Again, good agreement is observed and the relative error is less than 2%.

Table 2 :
Approximation errors in the derivatives of the function (, ) = sin  sin .
0 ,   =   −  0 , and   is weighing function with compact support.The solution of the derivatives is obtained by minimizing the norm , that is,