FINITE ELEMENT SOLUTION OF WAVE EQUATION WITH REDUCTION OF THE VELOCITY DISPERSION AND SPURIOUS REFLECTION

The objective of this paper is to present a finite element solution for the wave propagation problems with a reduction of the velocity dispersion and spurious reflection. To this end, a high-order two-step direct integration algorithm for the wave equation is adopted. The suggested algorithm is formulated in terms of two Hermitian finite difference operators with a sixth-order local truncation error in time. The two-node linear finite element presenting the fourth-order of local truncation error is considered. The numerical results reveal that although the algorithm competes with higher-order algorithms presented in the literature, the computational effort required is similar to the effort required by the average acceleration Newmark method. More than that, the integration with the lumped mass model shows similar results to the integration using the average acceleration Newmark for the consistent mass model, which involves a higher number of computational operations.


INTRODUCTION
In the numerical computation of wave equations, numerical wave velocity dispersion and spurious reflections for non-uniform meshes are a persistent problem arising from the inadequate discretization of the continuous. Wave propagation in discrete systems has been studied since the 17th century as described by Brillouin [1]. A one-dimensional lattice of point mass interconnected by springs has served, for example, as a model for sound wave propagation and propagation of waves in crystals. A lot of research effort has been devoted to dispersion and spurious reflections in the numerical integration of wave equations by finite element method Belytschko and Mullen [2], Liu, Sharan, and Yau [3] and Bazant [4]. Idesman [5] presented the optimal reduction in numerical dispersion for wave propagation problems using two-dimensional isogeometric elements.
Direct integration methods have been developed for a long time. One of the pioneers was the Houbolt method [6] using backward finite difference operators presenting numerical asymptotic annihilation and stability followed by the Newmark bi-parametric method [7] with numerical damping control. The single parametric Wilson method [8] was developed with improvements in the numerical damping but presenting initial displacement overshoot. Hilbert, Hughes, and Taylor [9] presented a new and efficient method with improvements in numerical dissipation. Bazzi and Anderheggen [10] developed a direct integration algorithm with new improved numerical dissipation and Hoff and Pahl [11] developed an improved variant of the Wilson method. Johnson [12] introduced the discontinuous Galerkin method. Finally, recently Soares Jr. [13] proposes a locally stabilized central difference method of four order accurate. The adopted high-order algorithm is formulated in terms of two Hermitian finite difference operators Collatz [14] with a sixth-order local truncation error. Because the developed algorithm considers the repeated differentiation of the governing equations, additional terms are required. Although the presence of these additional terms increases the number of the computational operation, the reduction obtained in the matrix factorization and higher orders of the relative radii errors is interesting attributes of the algorithm Laier [15]. More importantly, it should be emphasized that the results obtained with the lumped mass model for example are similar to those obtained by the Newmark method with the consistent mass model that requires greater computational effort.

ONE-DIMENSIONAL WAVE EQUATION
The classical one-dimensional wave equation can be written as (D'Alembert wave equation)

− = 
2 II c u u 0 (1) in which = ρ E c (2) and where c is the wave velocity of propagation, E is the modulus of elasticity, and ρ is the mass density. Roman numeral as exponent indicates space derivative and time derivative notation using upper dots are adopted.
Two main techniques are generally used to solve wave equation (1) (Clough, 1975). The first technique, referred to as integration in finite terms, is applied using normal mode superposition (a vibration solution), and the second technique, referred to as integration in non-finite terms, is known as the D´Alembert wave solution Achenbach [16].
The integration of equation (1) in non-finite terms (wave solution) is given using complex notation by where A and B are the amplitudes of the waves propagating in the positive and negative directions, respectively; β is the wave number, ω is the wave circular frequency, and i is the imaginary unit.

FINITE ELEMENT FORMULATION
The two-node element formulation, as illustrated in Figure 1, can be summarized by the following matrix equation Laier [17]: where  is the element length,  Note that equation (5) is a differential-difference equation (differential in time and difference in space) in which the space variable x is replaced by the discrete variable j where j=0,1… and  is the space increment.

LOCAL TRUNCATION ERROR OF THE NUMERICAL EQUILIBRIUM
Taking into account equation (1) and Taylor expansion of the function involved, the equation (5) Equation (6) indicates that both for the consistent mass model (

HIGH ORDER TIME INTEGRATION
To reduce wave velocity dispersion and spurious reflection in the numerical integration of the wave equation, the following coupled Hermitian high-order two-step version of the time integration algorithm suggested by Laier [15] is considered: where t ∆ is the time step. Equation (7) indicates that the coupled algorithm presents the sixth-order of convergence.
On the other hand, the two-step version of the average acceleration Newmark method presents the fourth-order of local truncation error.
Note that equation (7) is a difference equation involving four unknown functions ( u , u  , u  and u  ) in which the time variable t is replaced by the discrete variable k t ∆ where k=0,1… and t ∆ is the time increment.

VELOCITY DISPERSION ANALYSIS
The numerical integration presents two kinds of error: the first error called local error depends on the order of the approximation polynomial of the function to be integrated Smith [18], and the second called global error consists of the eigenvalue error and the eigenvector error. Although the solution of the wave equation (1) results in non-dispersive propagation (the velocity of propagation does not depend on the frequency), the numerical solution depends on the frequency of the wave resulting in dispersive propagation.
Equation (5) and its first-time derivative can be expressed as follows and equation (7) results in the following linear system of difference equations: Note that the equation (9) is a system of difference equations in space (discrete variable j ) and in time (discrete variable k t ∆ ) involving four unknown functions ( u , u  , u  and u  ), where where a is the time mesh parameter, b is the space mesh parameter, T is the wave period and cT λ = is the wavelength.
Note that the Courant number is the ratio b a (Belytshk0 and Mullen [2]).
The solution of the difference equation (9) for given wave circular frequency ω (a frequency domain analysis) can be expressed by Laier [17] ( By considering the equation (13) and substituting the equation (12) in the difference equation (9), the following eigenvalue problem is obtained after some algebraic manipulation: where i is the imaginary unit and After the mesh parameters, a and b are assumed (11), the following expression can be formulated: and the numerical wave number n β can be obtained from equation (14). The characteristic equation of the eigenvalue problem (14) in this case is expressed as where for lumped mass model ( It is important to emphasize that the numerical wave number n β  (the eigenvalue of equation (14)) for a given frequency is related to the exact wavenumber by where n c is the numerical wave velocity. Similarly, for a given wavenumber β , the eigenvalue of equation (14) is the numerical wave circular frequency n ω that is related to the exact wave circular frequency by n n c c ω = ω (20) Finally, it is important to mention that the eigenvalue problem (14) presents two excluding solutions. Table 1 presents the eigenvalue results (17) for relative numerical error for wave velocity for a given frequency and provides a comparison of the results for the relative numerical wave velocity error for the lumped mass model using the proposed high-order integration algorithm expressed by

WAVE DISPERSION RESULTS
with the results for the consistent mass model using the Newmark average acceleration method (results in brackets).
The numerical results indicate that for the fine mesh ( a 50 > and b 50 > ) the numerical wave velocity dispersion obtained by the adopted algorithm considering the lumped mass model are similar to the results obtained by the average acceleration Newmark method that considers the consistent mass model, for which the amount of operations is two times greater.

SPURIOUS WAVE REFLECTIONS
For a uniform mesh, the numerical integration of the wave equation (1) results in a dispersive wave, but with the same velocity of propagation throughout the entire domain for a given wave frequency. However, for irregular mesh, there is a reflection at the interface of elements with different lengths considering that there are different numerical impedances.
Consider a non-uniform finite element mesh in a configuration similar to the configuration shown in Figure 3. By considering an incident wave with unit amplitude traveling from left to right and arriving at the interface, spurious reflected wave and transmitted are generated.
where j=0,1,2…, tr A is the amplitude of the transmitted wave and tr β is the numerical wave number of the transmitted wave ( Figure 3). To calculate the spurious reflection amplitude re A and transmitted wave amplitude tr A , one should consider the displacement compatibility at the interface (node j), which is expressed by and the equilibrium at node j, which is written as Note that the following relationship between eigenvector components is also considered (see equation (14)) Note that re A is evaluated as a real component because the imaginary component disappears. Table 2 includes a comparison of results of the spurious wave reflection amplitude with the results of those by average acceleration Newmark method, considering the consistent mass model for an illustrative range of time and space mesh parameters, as presented in the brackets.
The results listed in Table 2 indicate that the magnitude of the relative error of the reflected spurious wave amplitude is the same as that obtained by the integration of the consistent mass model using the average acceleration Newmark method.

CONCLUSION
The wave velocity dispersion and spurious wave reflections at the interfaces of finite elements of different lengths (non-uniform mesh) predicted by the adopted high-order algorithm are examined here.
It should be emphasized that the numerical results indicate that the numerical wave velocity dispersion considering the lumped mass model is similar to the results of the integration by the average acceleration Newmark method for the fine mesh ( a 50 ≥ and b 50 ≥ ) considering the consistent mass matrix, for which the amount of computational operations is two times greater. In other words, the use of high-order direct integration algorithm can be highly advantageous in many cases. More than that, the magnitude of the spurious reflection also presents the same order of error as that obtained with the average acceleration Newmark integration for the consistent mass model.