An Enriched Finite Element Method with Appropriate Interpolation Cover Functions for Transient Wave Propagation Dynamic Problems

: A novel enriched ﬁnite element method (EFEM) was employed to analyze the transient wave propagation problems. In the present method, the traditional ﬁnite element approximation was enriched by employing the appropriate interpolation covers. We mathematically and numerically showed that the present EFEM possessed the important monotonic convergence property with the decrease of the used time steps for transient wave propagation problems when the unconditional stable Newmark time integration scheme was used for time integration. This attractive property markedly distinguishes the present EFEM from the traditional FEM for transient wave propagation problems. Two typical numerical examples were given to demonstrate the capabilities of the present method.

However, the conventional FE solutions for transient wave propagation problems always suffer from the numerical dispersion issue induced by the spatial discretization and temporal discretization [35][36][37][38][39]. The former is related to the spatial interpolation and the used mesh, and the latter is mainly dominated by the used temporal discretization scheme. As a result, the obtained numerical waves are always dispersive and inaccurate. Another important point is that the quality of the FE solutions for transient wave propagation analysis cannot be improved constantly by using the decreasing temporal discretization steps for a fixed mesh pattern [40][41][42][43]. In other words, the time step for temporal discretization should be carefully determined according to the used mesh and the considered wave speed. Note that many different wave components with different speeds are always involved in transient wave propagation problems; so in general, these different wave components cannot be solved simultaneously with good computation accuracy. To address the drawbacks of the standard FEM in solving the transient wave propagation problems, the relatively new numerical techniques have still received quite considerable research efforts in the computational mechanics field. This is also the main motivation of the present work.
In this paper, we mainly focus on employing a novel enriched finite element method (EFEM) to analyze the transient wave propagation problems. In the present EFEM, the constructed interpolation covers are used to enrich the original FE approximation to improve the performance of the standard FEM in transient wave propagation analysis. We show in detail that the monotonic convergence property with the decrescent temporal discretization steps can be achieved for transient wave propagation problems. In consequence, all different types of waves with different propagation speeds can be simulated accurately as long as the temporal discretization steps for time integration are small enough. We provide two typical numerical examples to examine and investigate the capacities of the present EFEM in a transient wave propagation analysis.
The main structure of the present paper is as follows: Section 2.1 gives the problem statement of this work. In Section 2.2, the formulation of the EFEM for the transient wave propagation problem is given. A complete dispersion analysis is performed for the EFEM in Section 2.3. In the following section, the performance of the EFEM in solving the transient wave propagation problem is evaluated by the typical numerical examples. Section 4 gives our conclusions.

Problem Statement
Consider a deformable elastic body occupying a bounded domain V. Assume that this elastic body is made of a homogeneous isotropic media and is subjected to a stress tensor τ(x, t) and body force vector f b (x, t). Based on the small displacement assumption, the strong from of this problem can be obtained as follows [1] ρ ..
in which u stands for the displacement, and ρ represents the mass density of the linear elastic media; .. u signifies two time derivatives.
By multiplying Equation (1) with a virtual displacement vector u and integrating by parts over the entire domain, we can obtain the following standard weak form of Equation (1) in which ε denotes the virtual strain, and f s f = τ · n is the imposed boundary traction vector. Using the standard finite element formulation, we can obtain the following matrix form of Equation (2) M ..
in which F(t) is the applied force vector, K is the stiffness matrix, M is the mass matrix and C is a matrix corresponding to the damping effects.

Formulation of the EFEM
In the enriched finite element model, the standard finite element interpolation for the considered field variable u at node i is enriched by the following expression [37,41]  in which a i1 a i2 a i3 · · · T denotes the additional unknown coefficients, and u i stands for the conventional nodal displacement corresponding to node i; L denotes the local interpolation basis functions. For the two-dimensional problems, which are discretized into regular mesh with average nodal space h, the following linear polynomials and trigonometric functions are always used as the local interpolation basis functions [35] in which the relative coordinate values(x i , y i ) are measured from node i, namely x i = x − x i and y i = y − y i ; the used trigonometric function is of q order; in this work, q = 1 and λ x = λ y = 2h are employed. Although the higher order of polynomials and trigonometric functions also can be used to construct the local enrichment functions, the involved numerical integration is accordingly more numerically expensive. Then Equation (4) can also be written by [37,41] From Equation (6), the distinct difference between the traditional FE approximation and the present enriched FE approximation can be seen because the present EFEM contains the additional interpolation cover in the interpolation scheme.
When the standard linear triangular element with the usual linear interpolation function N i is used, the global displacement u can be expressed by [37] in which j is the additional degree of freedom for node i, and m is the total number of nodes in the problem domain.
From Equations (6) and (7), it is seen that when we only use the constant term 1 in Equation (5) as the local enrichment function for required interpolation, the present EFEM will become the standard FEM. Therefore, the proposed EFEM in this work is actually the combination of the higher interpolation and the conventional FE interpolation. Using the present EFEM, the computation accuracy of the obtained numerical solutions can be markedly increased without adding the additional nodes, because the higher interpolation can be obtained using the low-order linear elements.

Dispersion Analysis
When the FEM-like numerical techniques are used for wave analysis, the involved numerical dispersion error issue should be carefully addressed because it is able to strongly influence the computation accuracy of the obtained numerical solutions. Therefore, the Mathematics 2022, 10, 1380 4 of 12 dispersion properties of the present EFEM will be evaluated carefully in this sub-section. As shown in Figure 1, in the dispersion analysis, the uniform triangular mesh is used, the blue line denotes the wave propagation direction θ and h stands for the average mesh size of the used mesh.

Dispersion Analysis
When the FEM-like numerical techniques are used for wave analysis, the involved numerical dispersion error issue should be carefully addressed because it is able to strongly influence the computation accuracy of the obtained numerical solutions. Therefore, the dispersion properties of the present EFEM will be evaluated carefully in this sub-section. As shown in Figure 1, in the dispersion analysis, the uniform triangular mesh is used, the blue line denotes the wave propagation direction θ and h stands for the average mesh size of the used mesh. By eliminating the time dependency from Equation (3) and without considering the boundary conditions, we have [35,37] ( ) in which [ ] lists all unknown solution coefficients, and k is the wave number.
Using the enriched finite element interpolation scheme shown in Equation (7) and referring to the uniform mesh shown in Figure 1, for one typical node, we can obtain [35,37] stiff mass is the amplitude vector for one typical node i [35,37]. By eliminating the time dependency from Equation (3) and without considering the boundary conditions, we have [35,37] in which a = u i a i1 a i2 a i3 · · · T lists all unknown solution coefficients, and k is the wave number.
Using the enriched finite element interpolation scheme shown in Equation (7) and referring to the uniform mesh shown in Figure 1, for one typical node, we can obtain [35,37] [D stiff − k 2 D mass ]a i = 0, (10) in which D stiff and D mass are the obtained Hermitian matrices of order n p × n p , corresponding to the stiffness matrix K and mass matrix M; n p is the total number of unknown solution coefficients for each node and a i = A 1 A 2 · · · A n p T is the amplitude vector for one typical node i [35,37]. If the nontrivial solutions to Equation (10) exist, we must have [35,37] det Note that the numerical wave number k h is involved in the matrices D stiff and D mass [35,37]; hence, the corresponding exact wave number k can be calculated using Equation (11) for any given k h . Owing to the numerical dispersion issue by the spatial discretization, k is in general different from k h . In this work, the obtained spatial dispersion error by the numerical methods is quantified by using the measure k/k h . Figure 2 shows the calculated spatial dispersion error k/k h from the standard FEM and the present EFEM as a function of the non-dimensional wave number kh along different wave propagation directions. It is seen that the spatial dispersion error from FEM is much larger than that from the present EFEM for the considered non-dimensional wave number. dispersion error by the numerical methods is quantified by using the measure h k k . Figure 2 shows the calculated spatial dispersion error h k k from the standard FEM and the present EFEM as a function of the non-dimensional wave number kh along different wave propagation directions. It is seen that the spatial dispersion error from FEM is much larger than that from the present EFEM for the considered non-dimensional wave number. The numerical dispersion error from the temporal discretization is also investigated in this section. Here, the additional boundary conditions and the damping effects are not The numerical dispersion error from the temporal discretization is also investigated in this section. Here, the additional boundary conditions and the damping effects are not considered. Assuming that all the variables in Equation (3) are time-harmonic, then the following matrix equation can be arrived at [35,37] M .. a + c 2 K a = 0, (12) in which c is the wave speed, a =âe j(k h n·x−ω h t) and ω h is the numerical angular frequency. Using the Newmark time integration method with the parameters δ = 1/2 and α = 1/4 for temporal discretization, we can obtain the following discretized equation [1] (13) in which ∆t is the used temporal discretization step, and ω is the exact angular frequency. Based on the formulation in Ref. [35], the following equation can be arrived at in which l = 1 + 1 4 ω 2 ∆t 2 , m = −2 + 1 2 ω 2 ∆t 2 and n = 1 + 1 4 ω 2 ∆t 2 . In Equation (14) we can see that the right side of the equation only contains the parameter ω∆t; hence, the parameter ω h ∆t is actually a function of the parameter ω∆t, namely where CFL = c∆t/h is the defined Courant-Friedrichs-Lewy (CFL) number [35,37]. Then, for transient wave propagations, the total numerical dispersion error can be calculated using the following equation Using the Taylor series expansion, we have Temporal dispersion error (17) Note , in which T and T h are, respectively, the exact and numerical period of the wave mode; then, we can have From above formulation, we can obtain the following important observations: (1) The total dispersion error c h /c for transient wave analysis clearly contains two different parts, namely k/k h and T/T h , which are respectively caused by the spatial discretization method (FEM and the present EFEM) and the time integration scheme for temporal discretization.
(2) T/T h is basically a monotonic function with respect to sufficiently small CFL numbers when the spatial discretization step h and the considered wave number k are determined. In other words, T/T h will approach 1 when the temporal discretization step ∆t, which is related to the CFL number, gets smaller, namely T/T h → 1 can be obtained when ∆t → 0 .
(3) Both the spatial and temporal discretization are able to affect the total dispersion error c h /c. The effects from these two different factors can be counteracted to some degree, because in general, T/T h 1 and k/k h 1. In addition, we also can observe that c h /c will tend to k/k h when ∆t → 0 (namely CFL → 0 ), because T/T h → 1 can be obtained when ∆t → 0 .
(4) There exists an optimal time step ∆t (namely the optimal CFL opt ) which can lead to almost non-dispersive numerical solutions, namely c h /c ≈ 1.
Based on the above findings, we can see that for transient wave analysis, the obtained total dispersion error c h /c from the proposed EFEM and the Newmark time integration scheme can be basically decreased by using the decreasing temporal discretization step ∆t.
As a result, the obtained numerical solutions from the EFEM and Newmark method will approach the "exact" solutions when the temporal discretization step ∆t tends to zero, because the involved spatial dispersion error k/k h is very small for kh π (see Figure 2); hence, the so-called monotonic convergence property with respect to the temporal discretization step ∆t can be reached, while this important and attractive property is not applicable for the standard FEM because the resultant spatial discretization error from the standard FEM is much larger than that from the EFEM (see Figure 2).
In addition, it should be noted that the involved linear dependence (LD) problem should be carefully addressed when the EFEM is employed for engineering computation. Owing to the LD problem, the nearly singular matrix equation can always be obtained, especially as the used local enrichment functions are constructed by the polynomials. In [37], the LD issue of the EFEM was investigated systematically, and the LD could be totally suppressed by an elegantly-designed scheme. Using this scheme, the relatively high computation accuracy of the EFEM could be completely maintained by removing the sufficiently few cover degrees of freedom (DOFs). In this work, the LD problem of the EFEM is addressed by using the scheme proposed in [37], and the performance of the EFEM for complicated transient waves is examined.

The Lamb's Problem
In this sub-section, the hardware configuration of the used laptop were two cores Intel 2.50 GHz CPU and 4 GB RAM. Firstly, the well-known Lamb's problem, with plane strain conditions in which the two-dimensional elastic wave propagations are involved, has been considered to investigate the capabilities of the present EFEM in transient wave propagation analysis. As shown in Figure 3, the material parameters of the considered elastic medium are mass density ρ =2200 kg/m 3 , P-wave velocity c p = 3200 m/s and S-wave velocity c s = 1848 m/s. The uniform triangular mesh with nodal space h = 20 m ias employed in the calculation. The following concentrated line load [37] is imposed at the free surface of the elastic medium.
where f p = 10 Hz is the peak frequency, and t s =0.1 s is the time shift.   Figure 4 gives the calculated displacement solutions from the standard FEM along the y-axis at the observation time t = 0.5 s; the higher and lower peaks correspond to the P-wave and S-wave, respectively. From the discussion in Section 2, it is known that the conventional FEM does not have the monotonic convergence property with the decreasing temporal discretization step; however, there exists an optimal time step (which is related to the crucial parameter CFL) for each wave component. Using the optimal time step, the minimal numerical dispersion error can be obtained [1]. In Figure  4, two different time steps Δ = t 0.004 s and Δ = t 0.007 s are used for time integration; these two time steps approximately correspond to the optimal CFL = 0.65 [40] for the P-wave and S-wave, respectively. From Figure 4, it is clearly seen that the P-wave and S-wave can be accurately simulated using the corresponding optimal steps; however, these two different elastic wave components with different speeds cannot be predicted simultaneously with sufficiently good accuracy.  Figure 4 gives the calculated displacement solutions from the standard FEM along the y-axis at the observation time t = 0.5 s; the higher and lower peaks correspond to the P-wave and S-wave, respectively. From the discussion in Section 2, it is known that the conventional FEM does not have the monotonic convergence property with the decreasing temporal discretization step; however, there exists an optimal time step (which is related to the crucial parameter CFL) for each wave component. Using the optimal time step, the minimal numerical dispersion error can be obtained [1]. In Figure 4, two different time steps ∆t = 0.004 s and ∆t = 0.007 s are used for time integration; these two time steps approximately correspond to the optimal CFL = 0.65 [40] for the P-wave and S-wave, respectively. From Figure 4, it is clearly seen that the P-wave and S-wave can be accurately Mathematics 2022, 10, 1380 8 of 12 simulated using the corresponding optimal steps; however, these two different elastic wave components with different speeds cannot be predicted simultaneously with sufficiently good accuracy.
conventional FEM does not have the monotonic convergence property with the decreasing temporal discretization step; however, there exists an optimal time step (which is related to the crucial parameter CFL) for each wave component. Using the optimal time step, the minimal numerical dispersion error can be obtained [1]. In Figure  4, two different time steps Δ = t 0.004 s and Δ = t 0.007 s are used for time integration; these two time steps approximately correspond to the optimal CFL = 0.65 [40] for the P-wave and S-wave, respectively. From Figure 4, it is clearly seen that the P-wave and S-wave can be accurately simulated using the corresponding optimal steps; however, these two different elastic wave components with different speeds cannot be predicted simultaneously with sufficiently good accuracy.  Figure 5 presents the corresponding numerical solutions from the present EFEM using three time steps for time integration. It is seen that both the P-wave and S-wave solutions will converge monotonically to the exact solutions with the decrease of the temporal discretization step. The reason for this is that the so-called monotonic convergence property can be basically reached when the transient wave propagations  Figure 5 presents the corresponding numerical solutions from the present EFEM using three time steps for time integration. It is seen that both the P-wave and S-wave solutions will converge monotonically to the exact solutions with the decrease of the temporal discretization step. The reason for this is that the so-called monotonic convergence property can be basically reached when the transient wave propagations are analyzed by the EFEM, so all considered wave components with different speeds can be accurately solved as long as the used time step is small enough. In addition, the von Mises stress distribution results from the present EFEM are also calculated and shown in Figure 6; we can see that the physical behaviors of all different wave components can be accurately simulated. From these numerical results, we can see that the proposed EFEM possesses very excellent properties and is a powerful numerical approach to tackle the transient wave propagations.

The Transient Wave Propagation along an Elastic Bar
In this sub-section, we further examine the abilities of the EFEM by studying the scalar wave propagation along a one-dimensional elastic bar. As shown in Figure 7, the right end of the bar is fixed, and the elastic bar is made of two different elastic media (which are represented in blue and black lines) with wave propagation speed c 1 = 1 m/s and c 2 = 2 m/s. The problem domain is discretized into a uniform mesh with nodal space h = 0.02 m. Assuming that a sinusoidal wave u = 0.8 sin(20πt)cm, t ∈ 0, 0.05 s begins to travel along this bar from the left end.
At the observation time t = 0.8 s, the displacement solutions from the standard FEM with linear elements and the present EFEM with linear polynomials and trigonometric enrichment functions are given in Figures 8 and 9. It is known that both reflected and transmitted waves will be induced in this problem, and it is easily to identify that in the figures, the small and large peaks, respectively, correspond to reflected and transmitted wave components. Similar to the previous numerical example, several different time steps are used in the calculation. From the results in Figures 8 and 9, the similar observations, which have been found in the previous numerical example, can again be obtained, namely all different wave components with different wave propagation speeds can be accurately predicted using the present EFEM when the used time step gets smaller, while the different waves cannot be predicted simultaneously in sufficiently good accuracy by the standard FEM.
Mathematics 2022, 10, x FOR PEER REVIEW 9 of 13 are analyzed by the EFEM, so all considered wave components with different speeds can be accurately solved as long as the used time step is small enough. In addition, the von Mises stress distribution results from the present EFEM are also calculated and shown in Figure 6; we can see that the physical behaviors of all different wave components can be accurately simulated. From these numerical results, we can see that the proposed EFEM possesses very excellent properties and is a powerful numerical approach to tackle the transient wave propagations.

The Transient Wave Propagation along an Elastic Bar
In this sub-section, we further examine the abilities of the EFEM by studying the scalar wave propagation along a one-dimensional elastic bar. As shown in Figure 7, the right end of the bar is fixed, and the elastic bar is made of two different elastic media are analyzed by the EFEM, so all considered wave components with different speeds can be accurately solved as long as the used time step is small enough. In addition, the von Mises stress distribution results from the present EFEM are also calculated and shown in Figure 6; we can see that the physical behaviors of all different wave components can be accurately simulated. From these numerical results, we can see that the proposed EFEM possesses very excellent properties and is a powerful numerical approach to tackle the transient wave propagations.

The Transient Wave Propagation along an Elastic Bar
In this sub-section, we further examine the abilities of the EFEM by studying the scalar wave propagation along a one-dimensional elastic bar. As shown in Figure 7, the right end of the bar is fixed, and the elastic bar is made of two different elastic media  At the observation time t = 0.8 s, the displacement solutions from the standard FEM with linear elements and the present EFEM with linear polynomials and trigonometric enrichment functions are given in Figures 8 and 9. It is known that both reflected and transmitted waves will be induced in this problem, and it is easily to identify that in the figures, the small and large peaks, respectively, correspond to reflected and transmitted wave components. Similar to the previous numerical example, several different time steps are used in the calculation. From the results in Figures 8 and 9, the similar observations, which have been found in the previous numerical example, can again be obtained, namely all different wave components with different wave propagation speeds can be accurately predicted using the present EFEM when the used time step gets smaller, while the different waves cannot be predicted simultaneously in sufficiently good accuracy by the standard FEM.

Conclusions
In this work, a novel enriched finite element method (EFEM) with additional appropriate interpolation cover functions was developed for transient wave propagation analysis. In the present EFEM, the additional linear polynomials and trigonometric functions were used to enrich the original nodal shape functions of the traditional FEM. It was mathematically and numerically shown that the present EFEM possessed the important and attractive monotonic convergence property with the decreasing temporal discretization steps in solving transient wave propagation problems; so, all different wave components with different propagation speeds could be accurately predicted using the sufficiently small time step. Therefore, the present EFEM is also a prospective powerful approach to deal with the elastic wave propagation in anisotropic media and multiple waves propagation in laminated composite structures.

Conclusions
In this work, a novel enriched finite element method (EFEM) with additional appropriate interpolation cover functions was developed for transient wave propagation analysis. In the present EFEM, the additional linear polynomials and trigonometric functions were used to enrich the original nodal shape functions of the traditional FEM. It was mathematically and numerically shown that the present EFEM possessed the important and attractive monotonic convergence property with the decreasing temporal discretization steps in solving transient wave propagation problems; so, all different wave components with different propagation speeds could be accurately predicted using the sufficiently small time step. Therefore, the present EFEM is also a prospective powerful approach to deal with the elastic wave propagation in anisotropic media and multiple waves propagation in laminated composite structures.

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