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.
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 around point remaining up to second-order terms yieldswhere , , , and .
By multiplying both sides of (1) with and ( 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 :where is a volume measure.
Repeating the same procedure with , , and instead of and , we get other three equations and the following is the result for : To approximate the integrations by Riemann sum, the volume of the support domain is divided into points with associated volumes , . Equations (2), (3), and the other three constitute a system of five equations written in matrix form aswithwhere with , , and , and 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 term (see (A.5)). If all the points in the domain have the same volume, at both sides of (4) will be cancelled, and the two final forms will be the same. However, 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 as in (1). Similar to the derivations of (2) and (3), but with different weight functions , , , , and , the following equations, as examples, for and , are obtained:Again the Riemann sum over the support domain is used to approximate the integrations and a system of five equations is obtained asCompared with formula (4), the only difference is the terms multiplied to both sides of (1). In GFDM, , , , , and are used instead of , , , , and 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: where can be , , , and 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 ( 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.
(a)
(b)
(a)
(b)
4.2. Two-Dimensional Case
For the functionits 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.
(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 arewhere is pressure, is particle velocity, is the density of the medium, and 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 and . That is, the derivatives of variable at point can be approximated by function values at points inside the support domain centered at aswhere and 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 haveBy taking both and in (12) as , 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 and 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) becomeswhere 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, , where 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 becomeswhere the particle velocity is a 2D vector .
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 and m/s and the time interval is set to be μ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 . 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 μs and 500 μs have good agreement with the FDTD solutions. Compared with FDTD, the relative errors concerning the pressure and are
(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. Figure 5 compares the solutions of FDTD and GFDTD at μs. The results along m are shown in Figure 5(c). Again, good agreement is observed and the relative error is less than 2%.
(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 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 m shown in Figure 6(c) indicates that, with irregular distribution of computational points, the Gaussian wave propagates as well as before.
(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 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(a). Due to the designed perturbation, the volume fluctuates around 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 The relatively small error is due to the gentle irregular point distribution. That is, if the point irregularity is small, the effect of is negligible. This is consistent with Benito’s results [9–11].
(a)
(b)
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 and are used as before and the time interval is μs. Inside the computational domain points with spatial interval mm are evenly spaced.
6.3.1. Source Term
The left is a wave source with pressure given by the Gaussian pulse:where and 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: where is the normal acoustic impedance on the boundary given byHere the normal sound absorption coefficient 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:By applying (12) to (23) and performing time integration, (23) degenerates toWhen GFDTD is used, the discrete form of (24) is obtained as
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.
(a)
(b)
(c)
(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.
Appendix
Conventional Derivation of GFDM
Considering the 2D case, for the same Taylor expansion in (1), we consider an energy norm :where , , , , and is weighing function with compact support.
The solution of the derivatives is obtained by minimizing the norm , that is, withFor example, the first equation is
Equation (A.4) and the other four give the following system:
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work is supported in part by the National Natural Science Foundation of China (nos. 51478305 and 61175016) and Key Program (no. 61233009). The authors thank the anonymous reviewers for their most useful suggestions.