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

Authors: Yingbin Chai; Yongou Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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


---

## Body

## 1. Introduction

The 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.Actually, a lot of different numerical techniques can be employed to deal with the wave problems, such as the finite element method (FEM) [1, 2], finite volume method (FVM) [3], the boundary element [4–6] or boundary-type numerical method [7–11], and meshless techniques [12–22].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–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–34]), including the smoothed finite element method (S-FEM) (see, e.g., [35–40]). The 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.The 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. The Bathe time integration scheme is a typical two-substep method. The 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. The 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. The 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.The 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 Section3; 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.


## 2. 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](1)u¨−c2∇2u=0,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):(2)MU¨+c2KU=0,in which U represents the unknown nodal variables; M and K are the obtained mass matrix and stiffness matrix, respectively:(3)M=∫ΩNTNdΩ,K=∫Ω∇NT∇NdΩ,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 edge-based 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.Figure 1
The description of forming the edge-based smoothing domains.In this work, the gradient smoothing operation is performed by smoothing the particle velocityv, and the smoothed particle velocity can be obtained by [49](4)v¯xk=∫ΩksvxkWx−xkdΩ,in which Ωks is the obtained smoothing domain for edge k and Wx−xk is a predefined smoothing function given by(5)Wx−xk=0,x∉Ωks,1Ask,x∈Ωks,in which Ask stands for the area of the smoothing domain.For the wave propagation in a prestressed membrane, the relationship between the particle velocityv and the displacement u is given by(6)v=−c∇u⋅n,in which n stands for the outward unit normal vector.Using (6) and (4), it can be rewritten by(7)v¯xk=−cAks∫Ωks∇udΩ=−cAks∫Γksu⋅ndΓ,in which Γks is the boundary of the corresponding smoothing domain.Following the conventional finite element interpolation steps, we have(8)v¯xk=−c∑i∈MkB¯ixui,in which Mk 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 ofNs segments; then smoothed gradient matrix B¯i can be obtained by(9)B¯ix=1Aks∑q=1Ns∑r=1NgwrNixnxny,in which Ng is the number of Gauss points in each segment and it is determined by the order of the used nodal shape functions, wr is the weighting coefficients, and nx and ny are the components of the outward unit normal vector.Once the smoothed gradient matrixB¯i is obtained, then the smoothed element stiffness matrix and smoothed global stiffness matrix can be obtained as in the standard FEM scheme:(10)K¯=∑K¯k=∑∫ΩksB¯TB¯dΩ=∑B¯TB¯Aks,in which K¯ is the smoothed global stiffness matrix and K¯k is the smoothed element stiffness matrix for edge k.


## 3. Dispersion Analysis

The 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 spaceh (see Figure 2).Figure 2
The employed uniform mesh pattern and the illustration of involved nodes for dispersion analysis. (a) FEM. (b) ES-FEM.
(a)(b)


### 3.1. 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]:(11)Ku−k2Mu=0,in which k is the wave number which is given by(12)k=ωc,in which ω is the angular frequency.For an interior nodeSp,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 Sp,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. The 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 Figure2 should have the same amplitude and can be represented by(13)u=Aejkhn⋅x,in which A is the amplitude of the numerical solution, j=−1, kh 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 Figure2 and substituting (13) into (11), we have(14)Dstiff−k2Dmassu=0,in which Dstiff and Dmass are the matrices corresponding to the system stiffness matrix K and system mass matrix M.For the standard FEM and ES-FEM models, the matricesDstiff and Dmass are obtained by(15)DstiffFEM=Kp,qFEM+Kp,q−1FEMe−jkhhcosθ+Kp,q+1FEMejkhhcosθ+Kp−1,q−1FEMejkhh−cosθ−sinθ+Kp−1,qFEMejkhh−sinθ+Kp+1,qFEMejkhhsinθ+Kp+1,q+1FEMejkhhcosθ+sinθ,(16)DmassFEM=Mp,qFEM+Mp,q−1FEMe−jkhhcosθ+Mp,q+1FEMejkhhcosθ+Mp−1,q−1FEMejkhh−cosθ−sinθ+Mp−1,qFEMejkhh−sinθ+Mp+1,qFEMejkhhsinθ+Mp+1,q+1FEMejkhhcosθ+sinθ,(17)DstiffES−FEM=Kp,qES−FEM+Kp,q−1ES−FEMe−jkhhcosθ+Kp,q+1ES−FEMejkhhcosθ+Kp−2,q−1ES−FEMejkhh−cosθ−2sinθ+Kp+2,q+1ES−FEMejkhhcosθ+2sinθ+Kp−1,q−2ES−FEMejkhh−2cosθ−sinθ+Kp−1,q−1ES−FEMejkhh−cosθ−sinθ+Kp−1,qES−FEMejkhh−sinθ+Kp−1,q+1ES−FEMejkhhcosθ−sinθ+Kp+1,q−1ES−FEMejkhh−cosθ+sinθ+Kp+1,qES−FEMejkhhsinθ+Kp+1,q+1ES−FEMejkhhcosθ+sinθ+Kp+1,q+2ES−FEMejkhh2cosθ+sinθ,(18)DmassES−FEM=Mp,qES−FEM+Mp,q−1ES−FEMe−jkhhcosθ+Mp,q+1ES−FEMejkhhcosθ+Mp−2,q−1ES−FEMejkhh−cosθ−2sinθ+Mp+2,q+1ES−FEMejkhhcosθ+2sinθ+Mp−1,q−2ES−FEMejkhh−2cosθ−sinθ+Mp−1,q−1ES−FEMejkhh−cosθ−sinθ+Mp−1,qES−FEMejkhh−sinθ+Mp−1,q+1ES−FEMejkhhcosθ−sinθ+Mp+1,q−1ES−FEMejkhh−cosθ+sinθ+Mp+1,qES−FEMejkhhsinθ+Mp+1,q+1ES−FEMejkhhcosθ+sinθ+Mp+1,q+2ES−FEMejkhh2cosθ+sinθ.To ensure that the nontrivial solutions to (14) exist, the following equation should be satisfied [46]:(19)detDstiff−k2Dmass=0.From (15)–(18), it is found that the matrices Dstiff and Dmass are the functions of the numerical wave number kh, so the relationship between the exact wave number k and the numerical wave number kh can be determined by (19). In other words, for any given numerical wave number kh, 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/kh.Figure3 gives the spatial discretization error k/kh in different wave propagation directions as a function of the normalized numerical wave number h/λh/2=khh/π (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. The 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. Therefore, 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.Figure 3
Comparisons of the spatial discretization error as a function of the normalized numerical wave number in different wave propagation directions for the standard FEM and the present ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)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 Figure3, 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.


### 3.2. Dispersion Error from Temporal Discretization

The 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:(20)u=Aejkn⋅x−ωht,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):(21)Dmassu¨+c2Dstiffu=0.Using the standard Bathe time integration scheme [45], (21) can be rewritten by(22)ut+2Δt+mut+Δt+nut=0,in which m=−288−94ω2Δt2/144+25ω2Δt2+ω4Δt4 and n=144+25ω2Δt2/144+25ω2Δt2+ω4Δt4, ω=kc is the exact angular frequency which can be obtained from (19) for a given numerical wave number kh, 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:(23)ωhΔt=arctan−4n−m2m,forΔt<Δt∗,arctan4n−m2m,forΔt>Δt∗,(24)ξh=−12lnnωhΔt,(25)AD=1−e−2πξh×100%,in which ξh is the numerical damping ratio, Δt∗ is the critical time step and it satisfies 4n−m2Δt=Δt∗=0, and AD denotes the percentage amplitude decay (AD) per period T.For any given numerical wave numberkh, 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(26)chc=ωh/khc=ωh/khω/k=kkhωhω=k/khTh/T,in which ch is the numerical wave speed and Th is the numerical period.From (26), we can find the following interesting points:(1)
The 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 byk/kh, is from the spatial discretization and the second part, which is defined by Th/T, is induced by the temporal discretization(2)
The spatial discretization errork/kh 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/kh will trend to 1)(3)
For the given wave speedc and nodal space h, the temporal discretization error Th/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, Th/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 ch/c) in various wave propagation directions as a function of the normalized wave number khh/π 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 Th/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. This 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.Figure 4
Comparisons of the total dispersion error results as a function of the normalized wave number in various wave propagation directions for both FEM and ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)Figure 5
Comparisons of the additional dispersion error results induced by the temporal discretization as a function of the normalized wave number in various wave propagation directions for both FEM and ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)Figure 6
Comparisons of the percentage amplitude decay (AD) as a function of the normalized wave number in various wave propagation directions for both FEM and ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)


## 3.1. 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]:(11)Ku−k2Mu=0,in which k is the wave number which is given by(12)k=ωc,in which ω is the angular frequency.For an interior nodeSp,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 Sp,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. The 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 Figure2 should have the same amplitude and can be represented by(13)u=Aejkhn⋅x,in which A is the amplitude of the numerical solution, j=−1, kh 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 Figure2 and substituting (13) into (11), we have(14)Dstiff−k2Dmassu=0,in which Dstiff and Dmass are the matrices corresponding to the system stiffness matrix K and system mass matrix M.For the standard FEM and ES-FEM models, the matricesDstiff and Dmass are obtained by(15)DstiffFEM=Kp,qFEM+Kp,q−1FEMe−jkhhcosθ+Kp,q+1FEMejkhhcosθ+Kp−1,q−1FEMejkhh−cosθ−sinθ+Kp−1,qFEMejkhh−sinθ+Kp+1,qFEMejkhhsinθ+Kp+1,q+1FEMejkhhcosθ+sinθ,(16)DmassFEM=Mp,qFEM+Mp,q−1FEMe−jkhhcosθ+Mp,q+1FEMejkhhcosθ+Mp−1,q−1FEMejkhh−cosθ−sinθ+Mp−1,qFEMejkhh−sinθ+Mp+1,qFEMejkhhsinθ+Mp+1,q+1FEMejkhhcosθ+sinθ,(17)DstiffES−FEM=Kp,qES−FEM+Kp,q−1ES−FEMe−jkhhcosθ+Kp,q+1ES−FEMejkhhcosθ+Kp−2,q−1ES−FEMejkhh−cosθ−2sinθ+Kp+2,q+1ES−FEMejkhhcosθ+2sinθ+Kp−1,q−2ES−FEMejkhh−2cosθ−sinθ+Kp−1,q−1ES−FEMejkhh−cosθ−sinθ+Kp−1,qES−FEMejkhh−sinθ+Kp−1,q+1ES−FEMejkhhcosθ−sinθ+Kp+1,q−1ES−FEMejkhh−cosθ+sinθ+Kp+1,qES−FEMejkhhsinθ+Kp+1,q+1ES−FEMejkhhcosθ+sinθ+Kp+1,q+2ES−FEMejkhh2cosθ+sinθ,(18)DmassES−FEM=Mp,qES−FEM+Mp,q−1ES−FEMe−jkhhcosθ+Mp,q+1ES−FEMejkhhcosθ+Mp−2,q−1ES−FEMejkhh−cosθ−2sinθ+Mp+2,q+1ES−FEMejkhhcosθ+2sinθ+Mp−1,q−2ES−FEMejkhh−2cosθ−sinθ+Mp−1,q−1ES−FEMejkhh−cosθ−sinθ+Mp−1,qES−FEMejkhh−sinθ+Mp−1,q+1ES−FEMejkhhcosθ−sinθ+Mp+1,q−1ES−FEMejkhh−cosθ+sinθ+Mp+1,qES−FEMejkhhsinθ+Mp+1,q+1ES−FEMejkhhcosθ+sinθ+Mp+1,q+2ES−FEMejkhh2cosθ+sinθ.To ensure that the nontrivial solutions to (14) exist, the following equation should be satisfied [46]:(19)detDstiff−k2Dmass=0.From (15)–(18), it is found that the matrices Dstiff and Dmass are the functions of the numerical wave number kh, so the relationship between the exact wave number k and the numerical wave number kh can be determined by (19). In other words, for any given numerical wave number kh, 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/kh.Figure3 gives the spatial discretization error k/kh in different wave propagation directions as a function of the normalized numerical wave number h/λh/2=khh/π (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. The 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. Therefore, 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.Figure 3
Comparisons of the spatial discretization error as a function of the normalized numerical wave number in different wave propagation directions for the standard FEM and the present ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)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 Figure3, 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.


## 3.2. Dispersion Error from Temporal Discretization

The 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:(20)u=Aejkn⋅x−ωht,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):(21)Dmassu¨+c2Dstiffu=0.Using the standard Bathe time integration scheme [45], (21) can be rewritten by(22)ut+2Δt+mut+Δt+nut=0,in which m=−288−94ω2Δt2/144+25ω2Δt2+ω4Δt4 and n=144+25ω2Δt2/144+25ω2Δt2+ω4Δt4, ω=kc is the exact angular frequency which can be obtained from (19) for a given numerical wave number kh, 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:(23)ωhΔt=arctan−4n−m2m,forΔt<Δt∗,arctan4n−m2m,forΔt>Δt∗,(24)ξh=−12lnnωhΔt,(25)AD=1−e−2πξh×100%,in which ξh is the numerical damping ratio, Δt∗ is the critical time step and it satisfies 4n−m2Δt=Δt∗=0, and AD denotes the percentage amplitude decay (AD) per period T.For any given numerical wave numberkh, 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(26)chc=ωh/khc=ωh/khω/k=kkhωhω=k/khTh/T,in which ch is the numerical wave speed and Th is the numerical period.From (26), we can find the following interesting points:(1)
The 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 byk/kh, is from the spatial discretization and the second part, which is defined by Th/T, is induced by the temporal discretization(2)
The spatial discretization errork/kh 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/kh will trend to 1)(3)
For the given wave speedc and nodal space h, the temporal discretization error Th/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, Th/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 ch/c) in various wave propagation directions as a function of the normalized wave number khh/π 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 Th/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. This 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.Figure 4
Comparisons of the total dispersion error results as a function of the normalized wave number in various wave propagation directions for both FEM and ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)Figure 5
Comparisons of the additional dispersion error results induced by the temporal discretization as a function of the normalized wave number in various wave propagation directions for both FEM and ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)Figure 6
Comparisons of the percentage amplitude decay (AD) as a function of the normalized wave number in various wave propagation directions for both FEM and ES-FEM. (a) FEM. (b) ES-FEM.
(a)(b)


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


### 4.1. Two-Dimensional Scalar Wave Propagation in a Square Domain

The first considered numerical example is the two-dimensional scalar wave propagation in a square prestressed membrane. This problem is described in Figure7. The 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. The computation domain is discretized into 2×80×80 uniform triangular elements with average nodal space h = 0.015 m and the time step Δt = 0.1 s is used for temporal discretization. For this scalar wave propagation problem, the transverse displacement u is governed by the following equation:(27)∇2u+Fcx=0,t=1c2∂2u∂t2.Figure 7
The problem description and used mesh patterns for the two-dimensional scalar wave propagation. (a) The geometry description of the square prestressed membrane. (b) The used uniform mesh.
(a)(b)The first concentrated force is a Ricker wavelet which is defined by [46](28)Fc=A1−2π2fp2t−ts2exp−π2fp2t−ts2,in which the magnitude A = 0.4 N, the peak frequency fp = 5 Hz, and the time shift ts = 0.25 s are used here.Figure8 gives the calculated displacement results from the standard FEM and the present ES-FEM along the wave propagation directions θ=45° at the observation time t = 0.9 s. The exact solution shown in the figure is obtained by using the Green function method. From the results shown in the figure, it is very clear to see that the present ES-FEM can produce better numerical solutions than the standard FEM and there are almost no spurious oscillations in the ES-FEM results. However, the FEM results deviate from the exact solution very much, and several spurious peaks can also be found. These findings indicate that the edge-based gradient smoothing technique used in the present ES-FEM is indeed effective to suppress the dispersion error, and then the accuracy of the calculated numerical solutions can be significantly improved.Figure 8
The calculated displacement results from the standard FEM and the present ES-FEM along the wave propagation directionsθ=45° for the first concentrated force case.Figures9 and 10 show the calculated displacement results from the standard FEM and the present ES-FEM along three different wave propagation directions (θ=0°, θ=22.5°, and θ=45°) at the observation time t = 0.9 s. It is seen that the numerical results from the standard FEM are strongly dependent on the wave propagation directions; namely, the standard FEM exhibits clear “numerical anisotropy” property in solving transient wave problems, while the present ES-FEM can successfully restrain this property to some degree. From Figure 10, we can find that the calculated solutions from the present ES-FEM are almost the same in different wave propagation directions. This numerical example demonstrates that the present ES-FEM works very well with the Bathe time integration scheme and indeed has obvious advantages in solving transient wave propagation problems compared to the standard FEM.Figure 9
Comparisons of the calculated displacement results from the standard FEM along various propagation directions for the first concentrated force case.Figure 10
Comparisons of the calculated displacement results from the present ES-FEM along various propagation directions for the first concentrated force case.Next, we consider another concentrated force which is defined by [46](29)Fc=1.6×102t0.1−t,0≤t<0.1,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.Figure 11
The calculated displacement results from the standard FEM and the present ES-FEM along the wave propagation directionsθ=45° for the second concentrated force case.Figure 12
Comparisons of the calculated displacement results from the standard FEM along various propagation directions for the second concentrated force case.Figure 13
Comparisons of the calculated displacement results from the present ES-FEM along various propagation directions for the second concentrated force case.


### 4.2. Two-Dimensional Scalar Wave Propagation Scattering Problem

The 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. This problem is described in Figure14. The external excitation load at the center of the membrane is still a Ricker wavelet with magnitude A = 0.4 N, the peak frequency fp = 10 Hz, and the time shift ts = 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.Figure 14
The problem description and used mesh patterns for the two-dimensional scalar wave scattering problem. (a) The geometry description of the considered problem. (b) The used mesh pattern.
(a)(b)At the observation timet = 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.Figure 15
Comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation directionθ=0° for two-dimensional scalar wave propagation scattering problem.Then, the corresponding displacement results along the other two different wave propagation directions (θ=22.5° and θ=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.Figure 16
Comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation directionθ=22.5° for two-dimensional scalar wave propagation scattering problem.Figure 17
Comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation directionθ=45° for two-dimensional scalar wave propagation scattering problem.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 Figures18 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.Figure 18
Snapshots of the displacement distributions from the FEM for the two-dimensional scalar wave scattering problem at different observation times: (a)t = 0.6 s; (b) t = 0.7 s; (c) t = 0.8 s; (d) t = 0.9 s.
(a)(b)(c)(d)Figure 19
Snapshots of the displacement distributions from the ES-FEM for the two-dimensional scalar wave scattering problem at different observation times: (a)t = 0.6 s; (b) t = 0.7 s; (c) t = 0.8 s; (d) t = 0.9 s.
(a)(b)(c)(d)


## 4.1. Two-Dimensional Scalar Wave Propagation in a Square Domain

The first considered numerical example is the two-dimensional scalar wave propagation in a square prestressed membrane. This problem is described in Figure7. The 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. The computation domain is discretized into 2×80×80 uniform triangular elements with average nodal space h = 0.015 m and the time step Δt = 0.1 s is used for temporal discretization. For this scalar wave propagation problem, the transverse displacement u is governed by the following equation:(27)∇2u+Fcx=0,t=1c2∂2u∂t2.Figure 7
The problem description and used mesh patterns for the two-dimensional scalar wave propagation. (a) The geometry description of the square prestressed membrane. (b) The used uniform mesh.
(a)(b)The first concentrated force is a Ricker wavelet which is defined by [46](28)Fc=A1−2π2fp2t−ts2exp−π2fp2t−ts2,in which the magnitude A = 0.4 N, the peak frequency fp = 5 Hz, and the time shift ts = 0.25 s are used here.Figure8 gives the calculated displacement results from the standard FEM and the present ES-FEM along the wave propagation directions θ=45° at the observation time t = 0.9 s. The exact solution shown in the figure is obtained by using the Green function method. From the results shown in the figure, it is very clear to see that the present ES-FEM can produce better numerical solutions than the standard FEM and there are almost no spurious oscillations in the ES-FEM results. However, the FEM results deviate from the exact solution very much, and several spurious peaks can also be found. These findings indicate that the edge-based gradient smoothing technique used in the present ES-FEM is indeed effective to suppress the dispersion error, and then the accuracy of the calculated numerical solutions can be significantly improved.Figure 8
The calculated displacement results from the standard FEM and the present ES-FEM along the wave propagation directionsθ=45° for the first concentrated force case.Figures9 and 10 show the calculated displacement results from the standard FEM and the present ES-FEM along three different wave propagation directions (θ=0°, θ=22.5°, and θ=45°) at the observation time t = 0.9 s. It is seen that the numerical results from the standard FEM are strongly dependent on the wave propagation directions; namely, the standard FEM exhibits clear “numerical anisotropy” property in solving transient wave problems, while the present ES-FEM can successfully restrain this property to some degree. From Figure 10, we can find that the calculated solutions from the present ES-FEM are almost the same in different wave propagation directions. This numerical example demonstrates that the present ES-FEM works very well with the Bathe time integration scheme and indeed has obvious advantages in solving transient wave propagation problems compared to the standard FEM.Figure 9
Comparisons of the calculated displacement results from the standard FEM along various propagation directions for the first concentrated force case.Figure 10
Comparisons of the calculated displacement results from the present ES-FEM along various propagation directions for the first concentrated force case.Next, we consider another concentrated force which is defined by [46](29)Fc=1.6×102t0.1−t,0≤t<0.1,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.Figure 11
The calculated displacement results from the standard FEM and the present ES-FEM along the wave propagation directionsθ=45° for the second concentrated force case.Figure 12
Comparisons of the calculated displacement results from the standard FEM along various propagation directions for the second concentrated force case.Figure 13
Comparisons of the calculated displacement results from the present ES-FEM along various propagation directions for the second concentrated force case.


## 4.2. Two-Dimensional Scalar Wave Propagation Scattering Problem

The 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. This problem is described in Figure14. The external excitation load at the center of the membrane is still a Ricker wavelet with magnitude A = 0.4 N, the peak frequency fp = 10 Hz, and the time shift ts = 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.Figure 14
The problem description and used mesh patterns for the two-dimensional scalar wave scattering problem. (a) The geometry description of the considered problem. (b) The used mesh pattern.
(a)(b)At the observation timet = 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.Figure 15
Comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation directionθ=0° for two-dimensional scalar wave propagation scattering problem.Then, the corresponding displacement results along the other two different wave propagation directions (θ=22.5° and θ=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.Figure 16
Comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation directionθ=22.5° for two-dimensional scalar wave propagation scattering problem.Figure 17
Comparisons of the calculated displacement results from the standard FEM and the ES-FEM along the wave propagation directionθ=45° for two-dimensional scalar wave propagation scattering problem.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 Figures18 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.Figure 18
Snapshots of the displacement distributions from the FEM for the two-dimensional scalar wave scattering problem at different observation times: (a)t = 0.6 s; (b) t = 0.7 s; (c) t = 0.8 s; (d) t = 0.9 s.
(a)(b)(c)(d)Figure 19
Snapshots of the displacement distributions from the ES-FEM for the two-dimensional scalar wave scattering problem at different observation times: (a)t = 0.6 s; (b) t = 0.7 s; (c) t = 0.8 s; (d) t = 0.9 s.
(a)(b)(c)(d)


## 5. Conclusions

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



---

*Source: 7180489-2020-07-31.xml*