Transient Wave Propagation Dynamics with Edge-Based Smoothed Finite Element Method and Bathe Time Integration Technique

In this work, the edge-based smoothed finite element method (ES-FEM) is incorporated with the Bathe time integration scheme to solve the transient wave propagation problems. +e edge-based gradient smoothing technique (GST) can properly soften the “overly soft” system matrices from the standard finite element approach; then, the spatial numerical dispersion error of the calculated solutions for wave problems can be significantly reduced. To effectively solve the transient wave propagation problems, the Bathe time integration scheme is employed to perform the involved time integration. Due to the appropriate “numerical dissipation effects” from the Bathe time integration method, the spurious oscillations in the relatively large wave numbers (high frequencies) can be effectively suppressed; then, the temporal numerical dispersion error in the calculated solutions can also be notably controlled. A number of supporting numerical examples are considered to examine the capabilities of the present approach. +e numerical results show that ES-FEM works very well with the Bathe time integration technique, and much more numerical solutions can be reached for solving transient wave propagation problems compared to the standard FEM.


Introduction
e wave propagation problems are usually encountered in real engineering applications. For several simple wave problems (such as a single wave propagating in a one-dimensional space), the exact solutions can be obtained using the analytical approach. When it comes to the relatively complex wave propagation problems, we have to resort to the numerical methods.
Among them, the classical finite element method is one of the most widely used and well-developed numerical approaches for analyzing wave problems. Based on the finite element method, several powerful and versatile commercial software packages have also been successfully developed for engineering computation. However, the finite element solutions for wave problems usually suffer from the numerical dispersion error issue for large wave numbers (namely, high frequencies) [23][24][25][26][27]. As a result, the calculated numerical solutions of wave problems from FEM are only reliable for relatively small wave numbers. In the large wave number range, the spurious oscillations always arise in the numerical results, and very inaccurate results can be frequently obtained. Quite importantly, when the low-order linear elements (such as triangular or quadrilateral element) are used, this troublesome dispersion issue will be even more severe.
In order to address the numerical dispersion error issue. Much research effort has been made and a variety of advanced or modified finite element schemes have been proposed (see, e.g., [28][29][30][31][32][33][34]), including the smoothed finite element method (S-FEM) (see, e.g., [35][36][37][38][39][40]). e S-FEM is developed by combining the classical finite element concepts and the generalized gradient smoothing technique (GGST) which is frequently employed in the meshless methods. Since the GGST can provide proper "softening effects" to the discretized model, then the inherent "overly stiff" property of the FEM model can be relieved to some extent and then the stiffness of obtained system matrices is closer to the real system. In consequence, the numerical dispersion error in solving wave propagation problems can be significantly reduced.
e present work mainly focuses on tackling the transient wave propagation problems. Many direct time integration techniques can be employed to solve the transient wave propagation dynamics, such as the central difference method [41], Wilson's θ method [42], Houbolt's method [43], Newmark's method [44], and so on. However, these different time integration schemes usually show several undesirable properties. For example, the central difference method is only a conditional stable method for transient dynamic analysis; the well-known Newmark method is unconditional stable scheme, but several extra adjustable parameters should be carefully determined for desirable computation accuracy; the Wilson θ method and Houbolt method always lead to a relatively large numerical error in period elongation and amplitude decay [45].
In this work, the Bathe implicit time integration technique is combined with the edge-based FEM (ES-FEM) to solve the transient wave propagation problems. e Bathe time integration scheme is a typical two-substep method.
e Newmark trapezoidal rule is used in the first substep and the 3-point Euler backward method is employed in the second substep. In the Bathe method, the appropriate numerical dissipation effects are introduced to suppress the spurious wave modes in high frequency range; then spurious oscillations in the calculated numerical solutions can be effectively eliminated and the solution accuracy can be significantly improved.
e numerical examples in this paper show that the edge-based smoothed FEM (ES-FEM) works very well with the Bathe implicit time integration technique in solving transient wave propagation problems. e numerical dispersion error in the calculated numerical results can be significantly suppressed and much more accurate numerical solutions can be obtained compared to the conventional FEM. e rest of the present work is structured as follows: in the next section the formulation of the ES-FEM for wave problems is briefly retrospected. A comprehensive dispersion analysis of the present ES-FEM with Bathe time integration scheme for transient wave propagation problems is given in Section 3; both the spatial dispersion error and temporal dispersion error are investigated in detail. A number of supporting numerical experiments and the related conclusions are then summarized in the remaining sections.

Formulation of the Edge-Based Smoothed FEM for Wave Problems
For a typical scalar wave problem (such as the wave propagation in a prestressed membrane), the governing equation is given by [46] € in which u denotes the field variable (such as the displacement), c is the wave speed, ∇ 2 stands for the Laplace operator, and the overdot represents the derivative with respect to time.
Following the basic finite element discretizations [47,48], we can obtain the corresponding matrix form of (1): in which U represents the unknown nodal variables; M and K are the obtained mass matrix and stiffness matrix, respectively: in which N is the used nodal shape function matrix. In this work, the edge-based smoothed FEM [49], in which the standard FEM is combined with the edgebased gradient smoothing technique, is used to deal with the transient wave propagation problems. In the ES-FEM model, the standard triangular mesh is used, and each edge-based smoothing domain is formed by sequentially connecting the endpoints of each edge and the centers of the neighboring elements (as shown in Figure 1). For the interior edges, two different elements will participate in constructing the smoothing domain, while for the edges on the global boundary, the corresponding smoothing domain only involves one element. In consequence, each edge corresponds to each smoothing domain uniquely and the gradient smoothing technique will be operated over these obtained smoothing domains.
In this work, the gradient smoothing operation is performed by smoothing the particle velocity v, and the smoothed particle velocity can be obtained by [49] v in which Ω s k is the obtained smoothing domain for edge k and W(x − x k ) is a predefined smoothing function given by in which A k s stands for the area of the smoothing domain. For the wave propagation in a prestressed membrane, the relationship between the particle velocity v and the displacement u is given by 2 Mathematical Problems in Engineering v � − c∇u · n, (6) in which n stands for the outward unit normal vector. Using (6) and (4), it can be rewritten by in which Γ s k is the boundary of the corresponding smoothing domain.
Following the conventional finite element interpolation steps, we have in which M k denotes the number of involved nodes which participate in forming the smoothed gradient matrix B i . In this work, the Gauss integration scheme is used to perform the related numerical integration along the boundary of smoothing domain which consists of N s segments; then smoothed gradient matrix B i can be obtained by in which N g is the number of Gauss points in each segment and it is determined by the order of the used nodal shape functions, w r is the weighting coefficients, and n x and n y are the components of the outward unit normal vector. Once the smoothed gradient matrix B i is obtained, then the smoothed element stiffness matrix and smoothed global stiffness matrix can be obtained as in the standard FEM scheme: in which K is the smoothed global stiffness matrix and K (k) is the smoothed element stiffness matrix for edge k.

Dispersion Analysis
e numerical solutions of the transient wave propagation problems usually suffer from the dispersion error issue. Both the spatial discretization and temporal discretization are able to lead to dispersion error. In this section, both the spatial discretization error and temporal discretization error will be investigated in detail by using the uniform node distributions with average nodal space h (see Figure 2).

Dispersion Error from Spatial Discretization.
For the time-independent form of the wave equation, the following matrix equation can be obtained without considering the boundary condition [47,48]: in which k is the wave number which is given by in which ω is the angular frequency.
For an interior node S (p,q) (here p and q represent the row number and column number, resp.), Figure 2 presents all the nodes (the blue solid nodes) which have contributions in forming the system matrix equations associated with the central node S (p,q) (the red solid node) for both FEM and ES-FEM models. We can see that in the present ES-FEM model more nodes are used to construct the system matrix equation compared to those in the standard FEM model. e reason for this is that the edge-based gradient smoothing operations are employed in the ES-FEM model.
Since no additional boundary conditions are considered here, so the numerical solutions corresponding to the involved nodes in Figure 2 should have the same amplitude and can be represented by in which A is the amplitude of the numerical solution, , k h is the numerical wave number which corresponds to the exact wave number k shown in (11), n is a unit vector denoting the wave propagation direction, and x is the position vector of the considered point.
By referring to the uniform mesh pattern shown in Figure 2 and substituting (13) into (11), we have in which D stiff and D mass are the matrices corresponding to the system stiffness matrix K and system mass matrix M.
For the standard FEM and ES-FEM models, the matrices D stiff and D mass are obtained by To ensure that the nontrivial solutions to (14) exist, the following equation should be satisfied [46]: From (15)-(18), it is found that the matrices D stiff and D mass are the functions of the numerical wave number k h , so the relationship between the exact wave number k and the numerical wave number k h can be determined by (19). In other words, for any given numerical wave number k h , we can obtain the corresponding exact wave number k by using (19); then we can obtain the spatial discretization error which is defined by k/k h . Figure 3 gives the spatial discretization error k/k h in different wave propagation directions as a function of the normalized numerical wave number h/(λ h /2) � k h h/π (here λ h is the numerical wavelength) for both the standard FEM and the present ES-FEM. It is found that both the FEM and ES-FEM can produce very small spatial discretization error for the relatively small normalized wave numbers, and the spatial discretization error will become larger with the increase of the considered normalized wave number. However, it is clearly seen that the spatial discretization error resulting from the present ES-FEM is much smaller than that from the standard FEM in the whole considered normalized wave number range. e reason for this is that the edge-based gradient smoothing operations in the ES-FEM are indeed helpful to suppress the spatial discretization error in solving wave propagation problems. erefore, it is reasonable to expect that the present ES-FEM can produce more accurate numerical solutions for solving transient wave propagation problems than the standard FEM.
In addition, we also can find that the spatial discretization error results from the standard FEM are strongly affected by the wave propagation directions; namely, the standard FEM shows clear "numerical anisotropy" property in solving wave propagation problems even if the uniform mesh pattern is used. From the results shown in Figure 3, it is seen that the present ES-FEM also suffers from the "numerical anisotropy" issue; however, it has been significantly relieved compared to the standard FEM.

Dispersion Error from Temporal Discretization.
e temporal discretization is also an indispensable ingredient for solving the transient wave propagation problems. Many direct time integration schemes can be used to deal with the transient wave propagation dynamics. However, the temporal discretization always can result in additional temporal discretization error which can degrade the quality of the calculated numerical solutions. In this section, we mainly focus on investigating the additional dispersion effects from the temporal discretization. Since the Bathe time integration technique is a very effective approach in dynamic analysis and shows several evident advantages compared to other time integration schemes [45], in this work the Bathe time integration method is used for temporal discretization.
If the additional boundary conditions are not considered, the fundamental solution to the matrix equation shown in (2) has the following form: in which ω h denotes the numerical angular frequency. By using (20) and referring to the node distributions shown in Figure 2, the following equation can be obtained from (2): Using the standard Bathe time integration scheme [45], (21) can be rewritten by t+2Δt u + m t+Δt u + n t u � 0, (22) in which m � − ((288 − 94ω 2 Δt 2 )/(144 + 25ω 2 Δt 2 + ω 4 Δt 4 )) and n � (144 + 25ω 2 Δt 2 )/(144 + 25ω 2 Δt 2 + ω 4 Δt 4 ), ω � kc is the exact angular frequency which can be obtained from (19) for a given numerical wave number k h , and Δt is the used time step for time integration. Based on the standard Bathe time integration scheme [45], the following equations can be directly obtained: for Δt > Δt * , in which ξ h is the numerical damping ratio, Δt * is the critical time step and it satisfies (4n − m 2 )| Δt�Δt * � 0, and AD denotes the percentage amplitude decay (AD) per period T. For any given numerical wave number k h , we can obtain the corresponding numerical angular frequency by using (23); then the total dispersion error for transient wave propagation problems can be obtained by in which c h is the numerical wave speed and T h is the numerical period. From (26), we can find the following interesting points: (1) e total dispersion error of the numerical solutions for transient wave propagation problems can be split into two different parts: the first part, which is defined by k/k h , is from the spatial discretization and the second part, which is defined by T h /T, is induced by the temporal discretization (2) e spatial discretization error k/k h is determined by (19) and it is a function of the average nodal space h; when the used nodal space h trends to zero, the obtained spatial discretization error will also trend to zero (namely, k/k h will trend to 1) (3) For the given wave speed c and nodal space h, the temporal discretization error T h /T is mainly determined by the used time step Δt; when the used time step Δt trends to zero, the obtained temporal discretization error will also trend to zero (namely, T h /T will trend to 1). For a fixed time step Δt � 0.5, Figure 4 gives the calculated total dispersion error results (which is measured by c h /c) in various wave propagation directions as a function of the normalized wave number k h h/π for both FEM and ES-FEM. It is found that the magnitude of the dispersion error from the standard FEM is clearly much larger than that from the present ES-FEM. Both the standard FEM and the present ES-FEM suffers from the "numerical anisotropy" issue; namely, the dispersion error results are different in different considered wave propagations; however, it is also found that this issue has been successfully suppressed to a certain extent by the present ES-FEM. In addition, the additional dispersion error induced by the temporal discretization (which is measured by T h /T) and the percentage amplitude decay (AD) per period T for the above-mentioned two different numerical techniques are displayed in Figures 5 and 6. Likewise, the very similar findings can be obtained from the figures. is means that the present ES-FEM with the Bathe time integration technique possesses better predictive capabilities in tackling the transient wave propagation  problems than the standard FEM, and much more accurate numerical solutions can be reached.

Numerical Examples
In this section, a number of supporting numerical examples are given to examine the effectiveness of the present ES-FEM with Bathe time integration scheme in solving transient wave propagation problems. It should be pointed out that the nonreflecting boundary conditions [50] (such as the Dirichlet-to-Neumann map, the absorbing boundary condition, and the perfect matched layer) are not considered in this work because all the wave components do not reach the boundary of the considered problem domain for the considered simulation time.

Two-Dimensional Scalar Wave Propagation in a Square
Domain. e first considered numerical example is the twodimensional scalar wave propagation in a square prestressed membrane. is problem is described in Figure 7. e length of the square problem domain is L � 2.4 m and the wave propagation speed is c � 1 m/s. A concentrated force is prescribed at the center of the square prestressed membrane. Due to symmetry, only a quarter of this in the involved problem domain is modelled for computation. e computation domain is discretized into 2 × 80 × 80 uniform e first concentrated force is a Ricker wavelet which is defined by [46] in which the magnitude A � 0.4 N, the peak frequency f p � 5 Hz, and the time shift t s � 0.25 s are used here.  Next, we consider another concentrated force which is defined by [46] F c � 1.6 × 10 2 t(0.1 − t), 0 ≤ t < 0.1, For this case, the displacements from the standard FEM and the present ES-FEM along the wave propagation direction θ � 45°at observation time t � 0.8 s are firstly calculated and shown in Figure 11. Similarly, we also investigated how the numerical solutions from the two mentioned methods are affected by the wave propagation directions, and the related numerical results are plotted in Figures 12 and 13. We again find that the present ES-FEM shows weaker "numerical anisotropy" property than the standard FEM in solving transient wave propagation problems, and better numerical results can be obtained.

Two-Dimensional Scalar Wave Propagation Scattering
Problem. e second considered numerical example is still the two-dimensional scalar wave propagation in a square prestressed membrane. In this case, four same holes are located uniformly in the problem domain. is problem is described in Figure 14. e external excitation load at the center of the membrane is still a Ricker wavelet with magnitude A � 0.4 N, the peak frequency f p � 10 Hz, and the time shift t s � 0.1 s. It is clear that both reflected wave and transmitted wave will arise when the original incident wave reach the boundary of the holes. In this numerical example, the abilities of the two different methods (FEM and ES-FEM) in predicting the behaviors of these different wave components will be investigated in detail. Likewise, only a quarter of the problem is modelled due to symmetry and the computation problem domain is discretized using standard triangular mesh with average nodal space h � 0.125 m, and the time step Δt � 0.1 s is used for temporal discretization. Note that the corresponding exact solutions to this problem are not easy to obtain; here the numerical solutions from the traditional FEM with a very fine mesh are provided as the reference solutions.
At the observation time t � 0.9 s, the comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation direction θ � 0°are plotted in Figure 15. In this wave propagation direction, there only exists the original incident wave, and the ES-FEM results are very close to the reference solution, while the FEM results are much worse than the ES-FEM results and many spurious peaks can be found in the FEM results. en, the corresponding displacement results along the other two different wave propagation directions (θ � 22.5°a nd θ � 45°) at the observation time t � 1 s are presented in Figures 16 and 17. In these wave propagation directions, both the reflected waves and the transmitted waves can be seen, and it is found that the present ES-FEM is indeed superior to the standard FEM in predicting the behaviors of the reflected waves and the transmitted waves because the ES-FEM results are much closer to the reference solutions than the FEM results.
Furthermore, several snapshots of the displacement distributions in the global computation domain from the standard FEM and ES-FEM at several observation times are shown in Figures 18 and 19. From the results shown in the figures, we can see that the dispersion error in the ES-FEM solutions is much smaller than that in the standard FEM solutions; hence, much smoother displacement distribution results can be reached and all the physical behaviors of the different wave components can be accurately predicted.

Conclusions
is paper focuses on presenting an edge-based smoothed FEM (ES-FEM) with the Bathe time integration technique for solving two-dimensional transient wave propagation problems. e dispersion error properties of the present approach for solving transient wave problems are investigated detailedly. It is found that the present ES-FEM works quite well with the Bathe time integration technique, and the dispersion error can be significantly reduced compared to the standard FEM. Quite importantly, the troublesome "numerical anisotropy" property of the standard FEM for transient wave problems can be clearly suppressed by the present ES-FEM; hence, the obtained numerical results are almost the same in the different wave propagation directions. e numerical results show that the present ES-FEM indeed surpasses the standard FEM in solving transient wave problems and can provide much more accurate numerical results in predicting the behaviors of the waves; hence, it can be regarded as a very promising numerical approach for solving practical transient wave dynamic problems.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.