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

Laier, José Elias

doi:10.1590/1679-78256221

ABSTRACT

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.

KEYWORDS:
finite element method; high-order time integration; local error control; wave velocity dispersion; spurious reflections

Graphical Abstract

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] Brillouin, L.; Wave propagation in periodic structures, Dover Publications, Inc., 2th Ed., 1953]. 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] T. Belytshko and R. Mullen, On dispersive properties of finite element solutions. In Modern Problems in Elastic Wave Propagation (Edited by J. Miklowitz and J. D. Achembach), John Wiley, New York, 67-78, 1978.], Liu, Sharan, and Yau [³[3] J. B. Liu, S. K. Sharan and L. Yao, Wave motion and its dispersive properties in finite element model with distortional elements, Computers & Structures, (1994), 52(2), 205-214.] and Bazant [⁴[4] Bazant, Z. P. Spurious reflections of elastic waves in nonuniform meshes of constant and linear strain finite elements. Computers & Structures, 11 (4), 451-459 1982]. Idesman [⁵[5] A. Idesman and B. Dey, Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements, Comput. Meth. Appl. Mech. Engng. (2017), 321, 235-268.] 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] Houbolt J. C., A recurrence matrix solution for the dynamic response of elastic aircraft, J, Aeronaut. Sci. 17, 540-550 (1950)] using backward finite difference operators presenting numerical asymptotic annihilation and stability followed by the Newmark bi-parametric method [⁷[7] Newmark, N. M. A method of computation for structural dynamics. Proc. Am. Soc. Civ. Engs, 85 (3), 67-94 (1959)] with numerical damping control. The single parametric Wilson method [⁸[8] Wilson E. L. A computer program for the dynamic stress analysis of underground structures, SESM Report 68-1, Department of Civil Engineering, University of California, Berkeley (1968)] was developed with improvements in the numerical damping but presenting initial displacement overshoot. Hilbert, Hughes, and Taylor [⁹[9] Hilber, H. M., Hughes, T.J.R. and Taylor, R.L. Improvement numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engrg. Structural dynamics, 5 283-292 (1977)] presented a new and efficient method with improvements in numerical dissipation. Bazzi and Anderheggen [¹⁰[10] Bazzi, G. and Anderheggen, E. The ρ-family of algorithms for time-step integration with improved numerical dissipation, Earthquake Engrg Structural Dynamics 10 537-550 (1982)] developed a direct integration algorithm with new improved numerical dissipation and Hoff and Pahl [¹¹[11] Hoff, C. and Pahl, P.J. Practical performance of the θ-method and comparison with other dissipative algorithms in structural dynamics Comput. Appl. Mech. Engrg 67 87-110 (1988)] developed an improved variant of the Wilson method. Johnson [¹²[12] Johnson C. Discontinuous Galerkin in finite element methods for second order hyperbolic problems Computer Methods in Applied Mechanics and Engineering 107 117-129 (1993).] 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] Collatz L. The numerical treatment of differential equations, Springer-Verlag, 2nd Ed., 1966.] 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] Laier J. E. Spectral analysis of a high-order Hermitian algorithm for structural dynamics, Appl. Math. Modelling, 2011, 35 n. 2, p. 965-971.]. 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)

c^{2} u^{I I} - \ddot{u} = 0

(1)

in which

c = \sqrt{\frac{E}{ρ}}

(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] Achenbach J. D. Wave propagation in elastic solids, North Holland, Amsterdam, 1973.].

The integration of equation (1) in non-finite terms (wave solution) is given using complex notation by

u = A \exp [i (β x - ω t)] + B \exp [i (β x + ω t)]

(3)

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] Laier JE. Mass lumping, dispersive properties and bifurcation of Timoshenko’s flexural waves. Adv. Eng. Software, 2007, 33, p. 605-610.]:

\frac{E}{l} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] {\begin{matrix} u_{j} \\ u_{j + 1} \end{matrix}} + ρ l [\begin{matrix} M_{1} & M_{2} \\ M_{2} & M_{1} \end{matrix}] {\begin{matrix} {\ddot{u}}_{j} \\ {\ddot{u}}_{j + 1} \end{matrix}} - {\begin{matrix} - N_{j} \\ N_{j + 1} \end{matrix}} = {\begin{matrix} 0 \\ 0 \end{matrix}}

(4)

where $l$ is the element length, $M_{1} = 1 / 3$ and $M_{2} = 1 / 6$ for the consistent mass model, $M_{1} = 1 / 2$ and $M_{2} = 0$ for the classical lumped mass model, and $N_{j} = E u_{j}^{I}$ and $N_{j + 1} = - E u_{j + 1}^{I}$ are normal forces at nodes j and j+1, respectively, as indicated in Figure 1.

Figure 1
Two node Finite Element

By considering equation (4), the equilibrium of a generic node j, as shown in Figure 2, which consists of the finite element version of the wave equation (1), can be expressed by

(- u_{j - 1} + 2 u_{j} - u_{j + 1}) + \frac{ρ l^{2}}{E} (M_{2} {\ddot{u}}_{j - 1} + 2 M_{1} {\ddot{u}}_{j} + M_{2} {\ddot{u}}_{j + 1}) = 0

(5)

Figure 2
Finite Element Equilibruium

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 l$ where j=0,1… and $l$ 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) permits to write [¹⁷[17] Laier JE. Mass lumping, dispersive properties and bifurcation of Timoshenko’s flexural waves. Adv. Eng. Software, 2007, 33, p. 605-610.] [¹⁸[18] Smith G.D. Numerical solution of partial differential equations: Finite Difference Methods, Clarenton Press - Oxford, 2nd Ed, 1978.]:

(- u_{j - 1} + 2 u_{j} - u_{j + 1}) + l^{2} (M_{2} u_{j - 1}^{I I} + 2 M_{1} u_{j}^{I I} + M_{2} u_{j + 1}^{I I}) = 0 + u_{j}^{I V} (M_{2} - 1 / 12) l^{4} + ...

(6)

Equation (6) indicates that both for the consistent mass model ( $5 u_{j}^{I V} l^{4} / 12$ ) and for the lumped mass model ( $- u_{j}^{I V} l^{4} / 12$ ) the local truncation error is fourth-order.

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] Laier J. E. Spectral analysis of a high-order Hermitian algorithm for structural dynamics, Appl. Math. Modelling, 2011, 35 n. 2, p. 965-971.] is considered:

\begin{array}{l} 12 (u_{k - 1} - 2 u_{k} + u_{k + 1}) + 6 Δ t ({\dot{u}}_{k - 1} - {\dot{u}}_{k + 1}) + Δ t^{2} ({\ddot{u}}_{k - 1} - 2 {\ddot{u}}_{k} + {\ddot{u}}_{k + 1}) = 0 + \frac{d^{6} u_{k}}{d t^{6}} \frac{Δ t^{6}}{12} + ... \\ 12 ({\dot{u}}_{k - 1} - 2 {\dot{u}}_{k} + {\dot{u}}_{k + 1}) + 6 Δ t ({\ddot{u}}_{k - 1} - {\ddot{u}}_{k + 1}) + Δ t^{2} ({\overset{⃛}{u}}_{k - 1} - 2 {\overset{⃛}{u}}_{k} + {\overset{⃛}{u}}_{k + 1}) = 0 + \frac{d^{7} u_{k}}{d t^{7}} \frac{Δ t^{6}}{12} + ... \end{array}

(7)

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$ , $\dot{u}$ , $\ddot{u}$ and $\overset{⃛}{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] Smith G.D. Numerical solution of partial differential equations: Finite Difference Methods, Clarenton Press - Oxford, 2nd Ed, 1978.], 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

(- {\dot{u}}_{j - 1} + 2 {\dot{u}}_{j} - {\dot{u}}_{j + 1}) + \frac{ρ l^{2}}{E} (M_{2} {\overset{⃛}{u}}_{j - 1} + 2 M_{1} {\overset{⃛}{u}}_{j} + M_{2} {\overset{⃛}{u}}_{j + 1}) = 0

(8)

and equation (7) results in the following linear system of difference equations:

\begin{array}{l} 12 (u_{k - 1} - 2 u_{k} + u_{k + 1}) + 6 Δ t ({\dot{u}}_{k - 1} - {\dot{u}}_{k + 1}) + Δ t^{2} ({\ddot{u}}_{k - 1} - 2 {\ddot{u}}_{k} + {\ddot{u}}_{k + 1}) = 0 \\ 12 ({\dot{u}}_{k - 1} - 2 {\dot{u}}_{k} + {\dot{u}}_{k + 1}) + 6 Δ t ({\ddot{u}}_{k - 1} - {\ddot{u}}_{k + 1}) + Δ t^{2} ({\overset{⃛}{u}}_{k - 1} - 2 {\overset{⃛}{u}}_{k} + {\overset{⃛}{u}}_{k + 1}) = 0 \\ (- u_{j - 1} + 2 u_{j} - u_{j + 1}) + ψ Δ t^{2} (M_{2} {\ddot{u}}_{j - 1} + 2 M_{1} {\ddot{u}}_{j} + M_{2} {\ddot{u}}_{j + 1}) = 0 \\ (- {\dot{u}}_{j - 1} + 2 {\dot{u}}_{j} - {\dot{u}}_{j + 1}) + ψ Δ t^{2} (M_{2} {\overset{⃛}{u}}_{j - 1} + 2 M_{1} {\overset{⃛}{u}}_{j} + M_{2} {\overset{⃛}{u}}_{j + 1}) = 0 \end{array}

(9)

Note that the equation (9) is a system of difference equations in space (discrete variable $j l$ ) and in time (discrete variable $k Δ t$ ) involving four unknown functions ( $u$ , $\dot{u}$ , $\ddot{u}$ and $\overset{⃛}{u}$ ), where

Ψ = \frac{ρ l^{2}}{E Δ t^{2}} = {(\frac{a}{b})}^{2}

(10)

with

\begin{array}{l} a = \frac{T}{Δ t} \\ b = \frac{λ}{l} \end{array}

(11)

where a is the time mesh parameter, b is the space mesh parameter, T is the wave period and $λ = c T$ is the wavelength. Note that the Courant number is the ratio $b / a$ (Belytshk0 and Mullen [²[2] T. Belytshko and R. Mullen, On dispersive properties of finite element solutions. In Modern Problems in Elastic Wave Propagation (Edited by J. Miklowitz and J. D. Achembach), John Wiley, New York, 67-78, 1978.]).

The solution of the difference equation (9) for given wave circular frequency $ω$ (a frequency domain analysis) can be expressed by Laier [¹⁷[17] Laier JE. Mass lumping, dispersive properties and bifurcation of Timoshenko’s flexural waves. Adv. Eng. Software, 2007, 33, p. 605-610.]

\begin{array}{l} u = A \exp [- i (β_{n} j l - ω k Δ t)] \\ \dot{u} = B \exp [- i (β_{n} j l - ω k Δ t)] \\ \ddot{u} = C \exp [- i (β_{n} j l - ω k Δ t)] \\ \overset{⃛}{u} = D \exp [- i (β_{n} j l - ω k Δ t)] \end{array}

(12)

where A, B, C, and D are the amplitudes of the corresponding waves and $β_{n}$ is the numerical wave number. The equations in (12) establish the following relationships (Euler´s Formula):

\begin{array}{l} u_{k - 1} - 2 u_{k} + u_{k + 1} = A \exp (- i β_{n} j l) \exp (i k ω Δ t) [2 \cos (ω Δ t) - 2] \\ u_{j - 1} - 2 u_{j} + u_{j + 1} = A \exp (- i β_{n} j l) \exp (i k ω Δ t) [2 \cos (β_{n} l) - 2] \\ u_{k - 1} - u_{k + 1} = A \exp (- i β_{n} j l) \exp (i k ω Δ t) [2 i \sin (ω Δ t)] \\ u_{j - 1} - u_{j + 1} = A \exp (- i β_{n} j l) \exp (i k ω Δ t) [2 i \sin (β_{n} l)] \end{array}

(13)

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:

\begin{array}{l} [\begin{matrix} 24 \cos (ω Δ t) - 24 & 12 i \sin (ω Δ t) & 2 \cos (ω Δ t) - 2 & 0 \\ 0 & 24 \cos (ω Δ t) - 24 & 12 i \sin (ω Δ t) & 2 \cos (ω Δ t) - 2 \\ 2 \cos (β_{n} l) - 2 & 0 & - 2 Ψ (M_{2} \cos (β_{n} l) + M_{1}) & 0 \\ 0 & 2 \cos (β_{n} l) - 2 & 0 & - 2 Ψ (M_{2} \cos (β_{n} l) + M_{1}) \end{matrix}] x \\ {\begin{matrix} A \\ B Δ t \\ C Δ t^{2} \\ D Δ t^{3} \end{matrix}} = {\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \end{matrix}} \end{array}

(14)

where i is the imaginary unit and $\cos (β_{n} l)$ is the eigenvalue. The corresponding eigenvector can be expressed as follows:

{\begin{matrix} 1 \\ \frac{(2 \cos (ω Δ t) - 2) [12 ψ (M_{2} \cos (β_{n} l) + M_{1}) + \cos (β_{n} l) - 1]}{12 ψ \sin (ω Δ t) (M_{2} \cos (β_{n} l) + M_{1})} i \\ \frac{\cos (β_{n} l) - 1}{ψ (M_{2} \cos (β_{n} l) + M {}_{1})} \\ \frac{(\cos (β_{n} l) - 1) (2 \cos (ω Δ t) - 2) [12 ψ (M_{2} \cos (β_{n} l) + M_{1}) + \cos (β_{n} l) - 1]}{12 ψ^{2} \sin (ω Δ t) {(M_{2} \cos (β_{n} l) + M_{1})}^{2}} i \end{matrix}}

(15)

After the mesh parameters, a and b are assumed (11), the following expression can be formulated:

\begin{array}{l} ψ = {(\frac{a}{b})}^{2} \\ ω Δ t = \frac{2 π}{a} \end{array}

(16)

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

d 1 \cos^{2} (β_{n} l) + d 2 \cos (β_{n} l) + d 3 = 0

(17)

where for lumped mass model ( $M_{1} = 1 / 2$ and $M_{2} = 0$ ) one has

\begin{array}{l} d 1 = 16 {(\cos (ω Δ t) - 1)}^{2} \\ d 2 = {(\cos (ω Δ t) - 1)}^{2} (96 ψ + 64) + 288 \sin^{2} (ω Δ t) \\ d 3 = {(\cos (ω Δ t) - 1)}^{2} [576 ψ^{2} - 96 ψ - 80] - 288 \sin^{2} (ω Δ t) \end{array}

(18)

It is important to emphasize that the numerical wave number $β_{n} l$ (the eigenvalue of equation (14)) for a given frequency is related to the exact wavenumber by

β_{n} = β \frac{c_{n}}{c}

(19)

where $c_{n}$ 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} = ω \frac{c_{n}}{c}

(20)

Finally, it is important to mention that the eigenvalue problem (14) presents two excluding solutions.

WAVE DISPERSION RESULTS

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

ε = \frac{c - c_{n}}{c}

(21)

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.

Thumbnail

Table 1
Relative numerical wave velocity error

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.

Figure 3
Non-uniform Finite Element Mesh

The incident wave be expressed as follows

u_{i n} = \exp [- i (ω k Δ t - β_{i n} j l)]

(22)

where j=0,-1,-2…, $β_{i n}$ is the numerical wave number of the incident wave (Jiang and Rogers 1991). Thus, the spurious reflected wave is given by

u_{r e} = A_{r e} \exp [- i (ω k Δ t + β_{i n} j l)]

(23)

where $A_{r e}$ is the amplitude of the spurious reflected wave. The transmitted wave is expressed as follows

u_{t r} = A_{t r} \exp [- i (ω k Δ t - β_{t r} j l)]

(24)

where j=0,1,2…, $A_{t r}$ is the amplitude of the transmitted wave and $β_{t r}$ is the numerical wave number of the transmitted wave (Figure 3).

To calculate the spurious reflection amplitude $A_{r e}$ and transmitted wave amplitude $A_{t r}$ , one should consider the displacement compatibility at the interface (node j), which is expressed by

A_{t r} = 1 + A_{r e}

(25)

and the equilibrium at node j, which is written as

\begin{array}{l} {\begin{matrix} - 1 & 1 \end{matrix}} {\begin{matrix} u_{- 1} \\ u_{0} \end{matrix}} + {\begin{matrix} \frac{1}{α} & - \frac{1}{α} \end{matrix}} {\begin{matrix} u_{0} \\ u_{+ 1} \end{matrix}} + ψ {\begin{matrix} M_{1} & M_{2} \end{matrix}} \frac{\cos (β_{i n} l) - 1}{M_{2} \cos (β_{i n} l) + M_{1}} {\begin{matrix} Δ t^{2} u_{- 1} \\ Δ t^{2} u_{0} \end{matrix}} + \\ ψ α {\begin{matrix} M_{1} & M_{2} \end{matrix}} \frac{\cos (β_{t r} L) - 1}{M_{2} \cos (β_{t r} L) + M_{1}} {\begin{matrix} Δ t^{2} u_{o} \\ Δ t^{2} u_{+ 1} \end{matrix}} = 0 \end{array}

(26)

in which

\begin{array}{l} α = \frac{L}{l} \\ u_{- 1} = [\exp (i β_{i n} l) + A \exp (- i β_{i n} l)] \exp (- i ω k Δ t) \\ u_{0} = (1 + A_{r e}) \exp (- i ω k Δ t) = A_{t r} \exp (- i ω k Δ t) \\ u_{+ 1} = A_{t r} \exp (- i β_{t r} L) \exp (- i ω k Δ t) = (1 + A_{r e}) \exp (- i β_{t r} L) \exp (- i ω k Δ t) \end{array}

(27)

Note that the following relationship between eigenvector components is also considered (see equation (14))

C Δ t^{2} = \frac{\cos (β l) - 1}{Ψ (M_{2} \cos (β l) + M_{1})} A

(28)

By substituting equation (25) in the equation (26) the following expression is obtained after some algebraic manipulation:

A_{r e} = \frac{(1 - P_{1}) \sin (β_{i n} l) + \frac{(1 - α P_{2})}{α} \sin (β_{t r} α l)}{(1 - P_{1}) \sin (β_{i n} l) + \frac{(α P_{2} - 1)}{α} \sin (β_{t r} α l)}

(29)

where

\begin{array}{l} P_{1} = \frac{M {}_{2}{(\cos (β_{i n} l) - 1)}}{M_{1} \cos (β_{i n} l) + M_{1}} \\ P_{1} = \frac{M {}_{2}{(\cos (β_{t r} α l) - 1)}}{M_{1} \cos (β {}_{t r}α l) + M_{1}} \end{array}

(30)

Note that $A_{r e}$ 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.

Thumbnail

Table 2
Relative spurious wave amplitude

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 \geq 50$ and $b \geq 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.

ACKNOWLEDGMENTS

The authors gratefully acknowledge the funding contributions to the Brazilian National Council for Technological and Scientific Development (CNPq).

REFERENCES

^[1]
Brillouin, L.; Wave propagation in periodic structures, Dover Publications, Inc., 2th Ed., 1953
^[2]
T. Belytshko and R. Mullen, On dispersive properties of finite element solutions. In Modern Problems in Elastic Wave Propagation (Edited by J. Miklowitz and J. D. Achembach), John Wiley, New York, 67-78, 1978.
^[3]
J. B. Liu, S. K. Sharan and L. Yao, Wave motion and its dispersive properties in finite element model with distortional elements, Computers & Structures, (1994), 52(2), 205-214.
^[4]
Bazant, Z. P. Spurious reflections of elastic waves in nonuniform meshes of constant and linear strain finite elements. Computers & Structures, 11 (4), 451-459 1982
^[5]
A. Idesman and B. Dey, Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements, Comput. Meth. Appl. Mech. Engng. (2017), 321, 235-268.
^[6]
Houbolt J. C., A recurrence matrix solution for the dynamic response of elastic aircraft, J, Aeronaut. Sci. 17, 540-550 (1950)
^[7]
Newmark, N. M. A method of computation for structural dynamics. Proc. Am. Soc. Civ. Engs, 85 (3), 67-94 (1959)
^[8]
Wilson E. L. A computer program for the dynamic stress analysis of underground structures, SESM Report 68-1, Department of Civil Engineering, University of California, Berkeley (1968)
^[9]
Hilber, H. M., Hughes, T.J.R. and Taylor, R.L. Improvement numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engrg. Structural dynamics, 5 283-292 (1977)
^[10]
Bazzi, G. and Anderheggen, E. The ρ-family of algorithms for time-step integration with improved numerical dissipation, Earthquake Engrg Structural Dynamics 10 537-550 (1982)
^[11]
Hoff, C. and Pahl, P.J. Practical performance of the θ-method and comparison with other dissipative algorithms in structural dynamics Comput. Appl. Mech. Engrg 67 87-110 (1988)
^[12]
Johnson C. Discontinuous Galerkin in finite element methods for second order hyperbolic problems Computer Methods in Applied Mechanics and Engineering 107 117-129 (1993).
^[13]
Soares Jr. D. A locally stabilized central difference method of fourth order accurate, Finite Element in Analysis and Design, 155 (2019) 1-10.
^[14]
Collatz L. The numerical treatment of differential equations, Springer-Verlag, 2^nd Ed., 1966.
^[15]
Laier J. E. Spectral analysis of a high-order Hermitian algorithm for structural dynamics, Appl. Math. Modelling, 2011, 35 n. 2, p. 965-971.
^[16]
Achenbach J. D. Wave propagation in elastic solids, North Holland, Amsterdam, 1973.
^[17]
Laier JE. Mass lumping, dispersive properties and bifurcation of Timoshenko’s flexural waves. Adv. Eng. Software, 2007, 33, p. 605-610.
^[18]
Smith G.D. Numerical solution of partial differential equations: Finite Difference Methods, Clarenton Press - Oxford, 2^nd Ed, 1978.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

Publication in this collection
11 Dec 2020
Date of issue
2020

History

Received
05 Aug 2020
Reviewed
08 Oct 2020
Accepted
14 Oct 2020
Accepted
15 Oct 2020

This is an open-access article distributed under the terms of the Creative Commons Attribution License

[1] ^[1]
Brillouin, L.; Wave propagation in periodic structures, Dover Publications, Inc., 2th Ed., 1953

[2] ^[2]
T. Belytshko and R. Mullen, On dispersive properties of finite element solutions. In Modern Problems in Elastic Wave Propagation (Edited by J. Miklowitz and J. D. Achembach), John Wiley, New York, 67-78, 1978.

[3] ^[3]
J. B. Liu, S. K. Sharan and L. Yao, Wave motion and its dispersive properties in finite element model with distortional elements, Computers & Structures, (1994), 52(2), 205-214.

[4] ^[4]
Bazant, Z. P. Spurious reflections of elastic waves in nonuniform meshes of constant and linear strain finite elements. Computers & Structures, 11 (4), 451-459 1982

[5] ^[5]
A. Idesman and B. Dey, Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements, Comput. Meth. Appl. Mech. Engng. (2017), 321, 235-268.

[6] ^[6]
Houbolt J. C., A recurrence matrix solution for the dynamic response of elastic aircraft, J, Aeronaut. Sci. 17, 540-550 (1950)

[7] ^[7]
Newmark, N. M. A method of computation for structural dynamics. Proc. Am. Soc. Civ. Engs, 85 (3), 67-94 (1959)

[8] ^[8]
Wilson E. L. A computer program for the dynamic stress analysis of underground structures, SESM Report 68-1, Department of Civil Engineering, University of California, Berkeley (1968)

[9] ^[9]
Hilber, H. M., Hughes, T.J.R. and Taylor, R.L. Improvement numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engrg. Structural dynamics, 5 283-292 (1977)

[10] ^[10]
Bazzi, G. and Anderheggen, E. The ρ-family of algorithms for time-step integration with improved numerical dissipation, Earthquake Engrg Structural Dynamics 10 537-550 (1982)

[11] ^[11]
Hoff, C. and Pahl, P.J. Practical performance of the θ-method and comparison with other dissipative algorithms in structural dynamics Comput. Appl. Mech. Engrg 67 87-110 (1988)

[12] ^[12]
Johnson C. Discontinuous Galerkin in finite element methods for second order hyperbolic problems Computer Methods in Applied Mechanics and Engineering 107 117-129 (1993).

[13] ^[13]
Soares Jr. D. A locally stabilized central difference method of fourth order accurate, Finite Element in Analysis and Design, 155 (2019) 1-10.

[14] ^[14]
Collatz L. The numerical treatment of differential equations, Springer-Verlag, 2^nd Ed., 1966.

[15] ^[15]
Laier J. E. Spectral analysis of a high-order Hermitian algorithm for structural dynamics, Appl. Math. Modelling, 2011, 35 n. 2, p. 965-971.

[16] ^[16]
Achenbach J. D. Wave propagation in elastic solids, North Holland, Amsterdam, 1973.

[17] ^[17]
Laier JE. Mass lumping, dispersive properties and bifurcation of Timoshenko’s flexural waves. Adv. Eng. Software, 2007, 33, p. 605-610.

[18] ^[18]
Smith G.D. Numerical solution of partial differential equations: Finite Difference Methods, Clarenton Press - Oxford, 2^nd Ed, 1978.

a	b
a	100	50	25	10
100	1.65 10-4	-6.59 10-4	-2.64 10-3	-1.69 10-2
100	(1.65 10-4)	(-3.28 10-4)	(-2.29 10-3)	(-1.57 10-2)
50	1.65 10-4	6.63 10-4	-2.64 10-3	-1.69 10-2
50	(1.15 10-3)	(6.58 10-4)	(-1.31 10-3)	(-1.47 10-2)
25	1.70 10-4	6.64 10-4	2.65 10-3	-1.69 10-2
25	(5.10 10-3)	(4.61 10-3)	(2.63 10-3)	(-1.08 10-2)
10	3.76 10-4	8.70 10-4	2.86 10-3	1.71 10-2
10	(3.29 10-2)	(3.24 10-2)	(3.04 10-2)	(1.66 10-2)

$α$	a=b=100	a=b=50	a=b=25	a=b=10
0.10	-2.444 10^-4	-9.791 10^-4	-3.940 10^-3	-2.574 10^-2
0.10	(-0.8149 10^-4)	(-3.268 10^-4)	(-1.320 10^-3)	(-0.8870 10^-2)
0.25	-2.314 10^-4	-9.272 10^-4	-3.732 10^-3	-2.444 10^-2
0.25	(-0.7717 10^-4)	(-3.094 10^-4)	(-1.250 10^-3)	(-0.8406 10^-2)
0.50	-1.852 10^-4	-7.420 10^-4	-2.990 10^-3	-1.974 10^-2
0.50	(-0.6174 10^-4)	(-2.476 10^-4)	(-1.001 10^-3)	(-0.6747 10^-2)
0.75	-1.082 10^-4	-4.331 10^-4	-1.749 10^-3	-1.170 10^-2
0.75	(-0.3602 10^-4)	(-1.445 10^-4)	(-0.5843 10^-3)	(-0.3958 10^-3)
1.10	5.136 10^-5	2.081 10^-4	8.379 10^-4	5.819 10^-3
1.10	(1.729 10^-5)	(0.6937 10^-4)	(2.809 10^-4)	(1.922 10^-3)
1.25	1.309 10^-4	5.580 10^-4	2.267 10^-3	1.590 10^-2
1.25	(0.4631 10^-4)	(2.859 10^-4)	(0.7532 10^-3)	(0.5182 10^-2)
1.50	3.091 10^-4	1.242 10^-3	5.065 10^-3	3.680 10^-2
1.50	(1.029 10^-4)	(0.4132 10^-3)	(1.677 10^-3)	(1.166 10^-2)
2.0	7.420 10^-4	2.990 10^-3	1.233 10^-2	9.929 10^-2
2.0	(2.471 10^-4)	(0.9928 10^-3)	(0.4044 10^-2)	(2.896 10^-2)

Brasil