Accurate finite-element modeling of wave propagation in composite and functionally graded materials

Composite Structures 117 (2014) 298–308

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Accurate finite-element modeling of wave propagation in composite and functionally graded materials A. Idesman ⇑ Department of Mechanical Engineering, Texas Tech University, Lubbock, TX 79409-1021, USA

a r t i c l e

i n f o

Article history: Available online 8 July 2014 Keywords: Waves Composite and functionally graded materials Time integration Finite elements Spurious oscillations

a b s t r a c t For the first time, we have obtained accurate numerical solutions for wave propagation in inhomogeneous materials under impact loading. We have extended the earlier developed numerical approach for elastodynamics problems in homogeneous materials to inhomogeneous materials. The approach includes the two-stage time-integration technique with the quantification and the filtering of spurious oscillations, the special design of non-uniform meshes as well as includes the standard finite elements and the elements with reduced dispersion. Similar to wave propagation in homogeneous materials in the 1-D case, we have obtained very accurate results for composite and functionally graded materials using the linear elements with lumped mass matrix and the explicit central difference method. We have also shown that specific non-uniform meshes yield much more accurate results compared to uniform meshes. We have also shown the efficiency of the finite elements with reduced dispersion compared with the standard finite elements. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction In the paper we will consider the accurate finite element modeling of wave propagation problems in inhomogeneous materials such as composite and functionally graded materials. The known finite element techniques for inhomogeneous materials (e.g., see [1–12]) are usually based on the explicit introduction of material properties as functions of the coordinates in the expressions for the mass and stiffness matrices. In some approaches, these properties are assumed to be constant within each element (but different for different elements), in other approaches they are approximated in terms of the standard finite element shape functions (the so-called graded elements). As shown in [4], the difference in numerical results between these formulations is not big. However, we have not seen in the literature the recommendations of how to select the dimensions of the finite elements depending on the variation of material properties. In many publications (e.g., see [1–12]), uniform meshes and meshes independent of material properties are used in calculations. However, the properties for composite and functionally graded materials can significantly differ in different locations, and an improper selection of the sizes of finite elements may lead to very inaccurate results or to a prohibitively large computation time for wave propagation problems ⇑ Tel.: +1 806 742 3563; fax: +1 806 742 3540. E-mail address: [email protected] http://dx.doi.org/10.1016/j.compstruct.2014.06.032 0263-8223/Ó 2014 Elsevier Ltd. All rights reserved.

if very fine meshes are used. Another issue with the finite element modeling of wave propagation in inhomogeneous material is related to the appearance and the quantification of spurious high-frequency oscillations. These oscillations can be very large (especially under impact loading) and can destroy the accuracy of numerical results. We should mention that accurate numerical results for wave propagation in inhomogeneous materials are very important and necessary for many engineering applications (e.g., see [13] for shock mitigation by the use of functionally graded materials). In our previous papers [14–21] we have developed the new numerical approach for wave propagation in homogeneous materials. This approach includes the two-stage time-integration technique with the stage of basic calculations as well as the filtering stage with the quantification and the filtering of spurious oscillations. This technique yields accurate numerical results for wave propagation problems at any loading (including impact loading) in the 1-D, 2-D and 3-D cases. The approach is implemented for implicit and explicit time-integration methods as well as for the low- and high-order standard finite elements, the spectral elements, isogeometric elements, and the linear elements with reduced dispersion. In this paper we will extend our numerical approach for elastodynamics problems to inhomogeneous materials. It will include the two-stage time-integration technique with the quantification and filtering of spurious oscillations, the special design of non-uniform meshes as well as will include the standard

A. Idesman / Composite Structures 117 (2014) 298–308

finite elements and the elements with reduced dispersion. Similar to wave propagation in homogeneous materials in the 1-D case, we will obtain very accurate results for composite and functionally graded materials using the linear elements with the lumped mass matrix and the explicit central difference method with the time increments equal to the stability limit. We will also show that specific non-uniform meshes yield much more accurate results compared to uniform meshes. For the first time, we will obtain accurate numerical solutions for wave propagation in inhomogeneous materials under impact loading. These numerical solutions do not include spurious oscillations and converge at mesh refinement (as shown below, existing approaches may yield divergent results at mesh refinement; see Figs. 12 and 14(a) and (b)). We will also show the efficiency of the linear finite elements with reduced dispersion compared with the standard linear finite elements. In our previous papers [15,20] we have shown that for elastodynamics problems the linear elements with reduced dispersion are more efficient (i.e., require less computation time at the same accuracy) than the standard quadratic elements (see [15]) and at moderate observation times they are more efficient than the high-order standard, spectral and isogeometric elements (see [20]). The paper consists of the extension of the two-stage timeintegration technique to wave propagation in inhomogeneous materials (Section 2.1), the design of non-uniform meshes for these problems (Section 2.2), and the examples of the accurate modeling of wave propagation problems in composite and functionally graded materials under impact loading (the most challenging case for accurate simulations) using the standard linear finite elements as well as the linear elements with reduced dispersion (Section 3). 2. Numerical technique The application of the space discretization to a system of partial differential equations for transient acoustics or transient linear elastodynamics leads to a system of ordinary differential equations in time

€ þ C U_ þ K U ¼ R; MU

ð1Þ

where M; C; K are the mass, damping, and stiffness matrices, respectively, U is the vector of the nodal displacement, R is the vector of the nodal load. Zero viscosity, C ¼ 0, is considered in the paper. Eq. (1) has the same form for homogeneous materials as well as for inhomogeneous materials including composites with a piecewise constant variation of material properties and functionally graded materials with a continuous variation of material properties. Due to the space discretization, the exact solution to Eq. (1) contains the numerical dispersion error; e.g., see [15,19,22–27] and many others. Therefore, even the exact time integration of Eq. (1) may lead to inaccurate results due to the space-discretization error. This can be seen for the problems with impact loadings for which large spurious oscillations may appear even for fine meshes in space; e.g., see [14–20]. In our papers for homogeneous materials [14–20] we have suggested several numerical techniques for the reduction of the numerical dispersion error of linear finite elements as well as for the accurate time integration of Eq. (1) without spurious oscillations under high-frequency and impact loadings. In this paper we will extend these techniques for wave propagation in inhomogeneous materials including composites and functionally graded materials. In the 1-D case, we will use the standard approach with linear 2-node finite elements for which the local consistent mass matrix M cons and the stiffness matrix K m for the mth finite element used in m Eq. (1) are given as (m ¼ 1; 2; 3; . . . ; n where n is the total number of finite elements):

M cons ¼A m

Z

dxm

0

qðsÞN T ðsÞNðsÞds

1 i2 h i s 1  dxsm qðsÞ  ðdxsm Þ2 qðsÞ dxm C ¼A Ads; i  2 0 s s 2 s  ð Þ q ðsÞ q ðsÞ dxm dxm dxm   Z dxm Z dxm 1 1 1 Km ¼ EðsÞBT ðsÞBðsÞds ¼ EðsÞds; 1 1 dxm 0 0 Z

dxm

299

0 h B @h

ð2Þ

ð3Þ

where N and B are the standard finite element shape and B matrices; A is the cross sectional area in the 1-D case; dxm is the length of the mth finite element. The local lumped mass matrix D can be obtained from the consistent mass matrix in Eq. (2) by the ‘‘row summation’’ technique. Along with the standard finite element formulations with the consistent and lumped mass matrices, we will use the finite element formulation with reduced dispersion based of a weighted average of the consistent M e and lumped D mass matrices with the weighting factor c (similar to that used in [15,19,25,26])

MðcÞ ¼ Dc þ M cons ð1  cÞ:

ð4Þ

This approach with reduced dispersion will be used below in the 1-D and 2-D cases (see Section 2.2). It is known that for homogeneous materials and uniform meshes, the weighting factor c ¼ 0:5 decreases the numerical dispersion error from the second order to the fourth order; see [15,19,25,26]. The value c ¼ 0:5 will be also used below in the numerical examples for inhomogeneous materials and non-uniform meshes. The global mass and stiffness matrices in Eq. (1) can be calculated by the standard summation of the corresponding local matrices. In Eqs. (2) and (3) we use the local Cartesian system with the origin at the left node. The difference between the finite element formulations for homogeneous and inhomogeneous materials consists in the fact that density q and Young’s modulus E in Eqs. (2) and (3) are constant for homogeneous materials and are functions of the space coordinate (the local coordinate s) for inhomogeneous materials. Therefore, in order to increase the accuracy of the results, some approaches were based on a variable density and Young’s modulus within a finite element where these quantities in Eqs. (2) and (3) were approximated with the help of the standard finite element shape functions; e.g., see [4,6,10,12]. However, numerical results showed that these modifications do not significantly change the accuracy of numerical results compared with a piecewise constant variation of density and Young’s modulus within the domain and the constant values of these parameters (calculated in the center of a finite element) within a finite element; see [4]. In our paper we will use the second possibility with a piecewise constant variation of density and Young’s modulus. This will allow us to extend some results for homogeneous materials to inhomogeneous materials (without the modifications of existing computer codes) and to obtain accurate numerical solutions for wave propagation in composite and functionally graded materials; see below. 2.1. Two-stage time-integration technique and its extension to composite and functionally graded materials The numerical solutions of wave propagation problems under high-frequency and impact loading include spurious high-frequency oscillations due to the large numerical dispersion error for high frequencies. In order to filter the spurious high-frequency oscillations, numerical dissipation (or artificial viscosity) is usually introduced for the time integration of Eq. (1). As we showed in our paper [14], the use of a time-integration method with numerical dissipation (or artificial viscosity) at each time increment leads to inaccurate numerical results for low frequencies as well, especially for a long-term integration. It is also unclear in this case

300

A. Idesman / Composite Structures 117 (2014) 298–308

how to select the amount numerical dissipation and how to quantify the range of spurious high frequencies to be filtered. To resolve these issues for wave propagation in homogeneous materials, we have developed the two-stage time-integration technique (see [14,17]) with the stage of basic computations and the filtering stage. This technique is based on the fact that for linear elastodynamics problems, there is no necessity to filter spurious oscillations at each time increment because the errors in high frequencies do not affect the accuracy of low frequencies during time integration; see [14]. In the current paper, at basic computations we use the standard explicit central-difference time-integration method (without numerical dissipation or artificial viscosity) or the standard implicit trapezoidal rule (without numerical dissipation or artificial viscosity) in order to obtain an accurate solution of the semi-discrete elastodynamics problem, Eq. (1) (this solution contains spurious high-frequency oscillations). For the filtering of spurious oscillations, the implicit time continuous Galerkin (TCG) method with large numerical dissipation developed in [14] is used at the filtering stage. For all elastodynamics problem, we use N ¼ 10 uniform time increments (5 positive plus 5 negative time increments) at the filtering stage. This means that there is no real time integration at the filtering stage (the sum of 10 time increments used at the filtering stage is zero). As shown in [14], this procedure is equivalent to the multiplication of each velocity and displacement of the uncoupled system of the semi-discrete  5 ð3þmÞ2 þX2 equations by a factor of ð3þmÞ (where X ¼ xj Dt and xj 2 þð2þmÞ2 X2 are the eigen-frequencies of the semi-discrete system, Dt is the time increment as well as m ¼ 15 is used) and does not require the modal decomposition and the calculation of eigen-frequencies. As can be seen, this factor is close to zero for large X and is close to unity for small X. The size Dt of time increments at the filtering stage indirectly defines the amount of numerical dissipation and the range of spurious oscillations and is calculated according to the following formulas (for uniform meshes and homogeneous materials)

  co T dx Dt ¼ a X0:1 ðNÞ; dx co

ð5Þ

with

   a2 co T co T ¼ a1 ; dx dx

a

in the 1-D case. Here, co ¼

ð6Þ qffiffiffi E

q

is the wave velocity; dx is the size of a

finite element in the 1-D case; T is the observation time; X0:1 ðN ¼ 10Þ ¼ 0:81 for the TCG method with N ¼ 10 time increments. Using the calibration procedure described in [17], the following coefficients a1 and a2 have been found for the linear finite elements: the standard approach with the consistent mass matrix – a1 ¼ 0:3574 and a2 ¼ 0:3204; the standard approach with the lumped mass matrix – a1 ¼ 0:3342 and a2 ¼ 0:3363; the approach with the averaged (c ¼ 0:5) mass matrix – a1 ¼ 0:3296 and a2 ¼ 0:2180. These coefficients a1 and a2 are calibrated in the 1-D case for the filtering of numerical results obtained at basic computations with very small time increments. We should also mention that the filtering stage can be applied in the beginning of calculations as a pre-processor, in the end of calculations as a post-processor or at some intermediate time (see [14–20] for numerous examples of the application of the two-stage time-integration technique). Eqs. (5) and (6) can be extended for composite materials as follows. Assume that a composite consists of G subdomains with different materials properties (the material properties within each subdomain are constant). Because each subdomain consists of a homogeneous material then a uniform mesh can be suggested for each subdomain for wave propagation problems (similar to wave

propagation in homogeneous materials; e.g., see [14–20]). As we also mentioned in [14], the range of spurious oscillations and the coefficients a1 and a2 are independent of the initial and boundary conditions. Therefore, in order to filter spurious high frequencies in each subdomain, the following time increments Dtp should be used for each subdomain at the filtering stage in the 1-D case:

 p  p c dx Dtp ¼ a op T X0:1 ðNÞ; cpo dx

ð7Þ

with



a

cpo pT dx



¼ a1



qffiffiffiffip

cpo pT dx

a2 ð8Þ

;

where cpo ¼ qE p and the superscript p corresponds to the pth subdomain of a composite (p ¼ 1; 2; . . . ; G). For wave propagation in composite isotropic materials in the 2-D and 3-D cases, compressional waves with the velocity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Ep ð1mp Þ cp1 ¼ qp ð1þ and shear waves with the velocity mp Þð12mp Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p Ep c2 ¼ 2qp ð1þmp Þ should be considered in each subdomain where mp is Poisson’s ratio in the pth subdomain (e.g., see our papers [14–20] for homogeneous materials). In this case, Eqs. (7) and (8) can be extended to the 2-D and 3-D problems as follows:

"

! p# cpm T dxj X0:1 ðNÞ p p m;j dxj cm " p #1a2 dxj a1 T a2 X0:1 ðNÞ ¼ max p m;j cm " # p 1a2 maxj dxj a1 T a2 X0:1 ðNÞ ¼ minm cpm  p 1a2 dxmax a1 T a2 X0:1 ðNÞ; ¼ max p cp2

Dtp ¼ max a

p ¼ 1; 2; . . . ; G;

ð9Þ

where cp2 ¼ minm cpm (m ¼ 1; 2) is the minimum value between the velocities of the compressional wave cp1 and the shear wave cp2 in p the pth subdomain (p ¼ 1; 2; . . . ; G), dxmax ¼ maxj dxj is the maximum dimension of finite elements along the axes xj in the pth subdomain (j ¼ 1; 2 for 2-D problems and j ¼ 1; 2; 3 for 3-D problems). At the filtering stage we currently use a time-integration method with the same time increments for the whole domain. Therefore, in order to filter the spurious oscillations in all subdomains, we should use the time increment equal to the maximum value among all Dt p ; i.e.,

Dt ¼ max Dtp ; p

p ¼ 1; 2; . . . ; G;

ð10Þ

where Dt p is calculated according to Eqs. (7) and (8) in the 1-D case and to Eq. (9) in the multidimensional case. We should mention that for the subdomains with Dt p < Dt, the range of the filtered high frequencies will be greater than necessary; i.e., some actual high frequencies will be filtered for these subdomains. However, the global numerical solution does not include spurious oscillations after the filtering stage and the numerical results converge to an exact solution at mesh refinement. Eqs. (7)–(10) can be easily extended to the wave propagation in functionally graded materials. As we mentioned above, in our approach we use the constant material properties within each finite element. In this case, we can use Eq. (10) with the numbers of subdomains G equal to the total number n of finite elements; i.e., G ¼ n. The material properties in each element are calculated in the center of the element.

301

A. Idesman / Composite Structures 117 (2014) 298–308

2.2. Finite element meshes for wave propagation in composite and functionally graded materials 2.2.1. Composite materials First we will discuss the selection of finite element meshes in the 1-D case. Waves in isotropic homogeneous materials propagate in all directions and in all locations with the same velocity co . Therefore, it seems that the selection of uniform finite element meshes for such problems is rational. For composite materials, the wave velocity cp0 is different for different subdomains and the selection of optimal finite element meshes is not obvious. Let us discuss some motivations for the selection of finite element meshes for composite materials. As we explained in Section 2.1, in each subdomain of a composite material, a uniform mesh can be recommended. In this case the range of spurious frequencies in a finite element solution for each subdomain can be indirectly estimated in terms of the time increments Dt p given by Eqs. (7) and (8). In order to filter only spurious frequencies in all subdomains at the filtering stage with the procedure described by Eqs. (7)–(10), the time increments Dtp given by Eqs. (7) and (8) should be the same for all subdomains (p ¼ 1; 2; . . . ; G); i.e.,

Dt p ¼ Dt ¼ const for p ¼ 1; 2; . . . ; G:

ð11Þ

In this case the range of actual frequencies is the same for all subdomains and no actual frequencies are filtered at the filtering stage; see the comments after Eq. (10). According to Eq. (7), we can rewrite Eq. (11) as follows:

cpo p ¼ const for dx

p ¼ 1; 2; . . . ; G;

ð12Þ

i.e., the dimensions of finite elements in each subdomain should be proportional to the wave velocity cpo in this subdomain. The same condition, Eq. (12), for the selection of finite element meshes for composite materials can be obtained if we use the explicit central-difference method for the time integration and a uniform mesh with linear finite elements for each subdomain. It is known that in the case of homogeneous materials, the integration of Eq. (1) by the explicit central difference method with the time increments equal to the stability limit Dt ¼ co dx yields exact solutions at the stage of basic computations and the filtering stage is not necessary. With the meshes that meet the condition given by Eq. (12), the stability limit is the same for all subdomains and the integration with the time increments equal to the stability limit yields exact solutions for composite materials at the stage of basic computations and the filtering stage is not necessary (see below the numerical examples). Eqs. (11) and (12) can be extended to the multidimensional case as follows. Eq. (11) is used without any modification. Then, according to Eq. (9), we will get the following condition for the selection of optimal finite element meshes in the multidimensional case

cp2 ¼ const for p dxmax

p ¼ 1; 2; . . . ; G;

ð13Þ

i.e., the maximum dimensions of finite elements in each subdomain should be proportional to the shear wave velocity cp2 in this subdomain. 2.2.2. Functionally graded materials (1-D case) The finite element simulation of wave propagation in functionally graded materials can be considered as a particular case of the simulation in a composite material. In this case, each finite element represents a subdomain with its own material properties; i.e., we can use Eq. (12) with G ¼ n where n is the total number of elements. Let us consider the calculation of an 1-D non-uniform mesh for functionally graded materials that yields (according to Eq. (12)) p co Dt the same Courant number dx p for all finite elements. The wave

velocity co ¼ f ðxÞ is a given function of f ðxÞ for a continuous problem and it is considered as a constant within an element cpo ¼ f ð xp Þ for a discretized problem where  xp is a given point within a finite element p (p ¼ 1; 2; . . . ; n). For the same Courant number, the p length lp ¼ dx of each element should be proportional to the wave p velocity co ; i.e.,

lp ¼ q cpo ;

ð14Þ

where q is a constant to be determined for a given total number n of finite elements. The wave velocity cpo is calculated at point xp within the pth finite element; i.e.,

x1 ¼ d1 l1

and xp ¼

p1 X lm þ d p lp ;

p ¼ 2; 3; . . . ; n;

ð15Þ

m¼1

where dp (0 6 dp 6 1) are the given coefficients. Combining Eqs. (14) and (15) we have the following system for the calculation of the coefficient q for a given number n of elements:

! p1 X l1 ¼ qf ðd1 l1 Þ and lp ¼ qf lm þ dp lp ;

p ¼ 2; 3; . . . ; n;

ð16Þ

m¼1

L¼

n X lm ;

ð17Þ

m¼1

where L is the total length of the bar. The non-linear system of Eqs. (16) and (17) can be iteratively solved for the determination of q and lp . However, in some cases the solution of Eqs. (16) and (17) can be found analytically. For example, for the linear variation of the wave velocity along the bar c ¼ f ðxÞ ¼ ax þ b (a and b are the given coefficients) and for dp ¼ 0:5 for p ¼ 1; 2; . . . ; n (i.e., the wave velocity in each element is calculated in the center of this element; see Eq. (15)) Eq. (16) reduces to " ! # p1 X ð18Þ l1 ¼ qð0:5al1 þ bÞ and li ¼ q a lm þ 0:5lp þ b ; p ¼ 2;3;...; n: m¼1

From Eq. (18) it follows that

ðlp  lp1 Þ



 1  0:5 ¼ lp1 ; aq

p ¼ 2; 3; . . . ; n

ð19Þ

or

  2aq lp1 ¼ rlp1 ; li ¼ 1 þ 2  aq

p ¼ 2; 3; . . . ; n;

ð20Þ

2aq with r ¼ 1 þ 2aq ¼ 2þaq ; i.e., lp is a geometric progression and 2aq

L¼

n X 1  rn lm ¼ l1 : 1r m¼1

ð21Þ

2qb Then, inserting l1 ¼ 2qa ¼ ðr  1Þ ba (see also Eq. (18)) into Eq. (21) we can find

r¼



aL þ1 b

1n ;

q¼

2ðr  1Þ ; aðr þ 1Þ

l1 ¼

2qb 2  aq

ð22Þ

for the given material parameters a and b, the given length L of the bar, and the given total number n of finite elements. Then, using r and l1 given by Eq. (22) we can find the lengths of all finite elements lp ¼ rlp1 (p ¼ 2; 3; . . . ; n) and the corresponding nodal coordinates. Remark. The generation of optimal finite element meshes for functionally graded materials in the multidimensional case is more difficult compared with the 1-D case. This will be considered elsewhere. As a possible solution, transition regions (similar to those used below in Section 3.3) can be used for the generation of optimal finite element meshes.

302

A. Idesman / Composite Structures 117 (2014) 298–308

3. Numerical examples In order to show the applicability of the new approach to composite materials it is sufficient to consider a composite made of two materials with different wave velocities. Below we will consider wave propagation under impact loading in a 1-D bi-material bar, in a 1-D bar made of functionally graded materials and in a 2-D bi-material cylinder. We will show that for the time integration of Eq. (1) in the 1-D case by the explicit central-difference method with the time increments equal to the stability limit and for special non-uniform meshes with standard linear finite elements, the numerical solutions to these problems are close to the exact solutions and do not include spurious oscillations. For other values of time increments, other explicit and implicit time-integration methods or other meshes, the numerical solutions to these problems include spurious high-frequency oscillations that increase with the increase in the observation time. Therefore, the two-stage time-integration technique described in Section 2.1 is used in such cases. The filtering stage of this technique indirectly quantifies and filters spurious high frequencies and is based on the TCG method with N ¼ 10 time increments (5 positive plus 5 negative time increments) the size of which is calculated according to Eqs. (7)–(9) for all examples presented below. For convenience, we will use the following dimensionless quantities for the numerical examples considered below:

L¼

Ld ; L0

x¼

E¼

Ed ; E0

q¼

xd ; L0

t¼

t d co ; L0

u¼

ud co ; L0 v 0

v¼

vd ; v0

qd q0

ð23Þ

for the 1-D problems and

Rd Ld rd zd td co ; L¼ ; r¼ ; z¼ ; t¼ ; R0 R0 R0 R0 R0 d ud co vd E qd ui ¼ i ; v i ¼ i ; E ¼ ; q ¼ R0 v 0 v0 E0 q0

R¼

ð24Þ

for the 2-D problems. Here, L (Ld ), x (xd ), t (td ), u (ud ), v (v d ) are the dimensionless (actual) length, coordinate, time, displacement and velocity in the 1-D case; L (Ld ), R (Rd ), r (rd ), z (zd ), t (t d ), ui (ud ), v (v di ) are the dimensionless (actual) length and radius of a cylinder, coordinates, time, displacements and velocities (i ¼ r; z) in the 2-D case; L0 is the given length, R0 is the given radius of a cylinder, v 0 is the given velocity applied at the left end of the bar for the problems considered below, E0 and q0 is the actual Young’s modulus and density ratio at the left end of the bar (at x ¼ 0) or cylinder (at qffiffiffiffi z ¼ 0), co ¼ qE0 . 0

3.1. Wave propagation in a 1-D bi-material bar under impact First, wave propagation in a 1-D bi-material bar of length L ¼ 4 under impact loading is considered (see Fig. 1). The bar is made of two materials with the following properties: Young’s modulus is chosen to be E ¼ 1 for 0 6 x 6 2 and E ¼ 16 for 2 6 x 6 4; the density to be q ¼ 1 for the whole bar where the origin of the Cartesian system is located at the left end of the bar. This corresponds to the following wave velocity: co ¼ 1 for 0 6 x 6 2 and co ¼ 4 for

Fig. 1. A bi-material bar under impact loading.

2 6 x 6 4. The following boundary conditions are applied: the displacement uð0; tÞ ¼ t (which corresponds to the velocity v ð0; tÞ ¼ v 0 ¼ 1) and uð4; tÞ ¼ 0 (which corresponds to the velocity v ð4; tÞ ¼ 0). Initial displacements and velocities are zero; i.e., uðx; 0Þ ¼ v ðx; 0Þ ¼ 0. The observation time is chosen to be T ¼ 98:3. During this long time the velocity pulse travels within the bar with multiple reflections from each end of the bar and from the interface between two materials at x ¼ 2. The analytical solution to this problem includes a piecewise constant distribution of the velocity along the bar at any time. The numerical solutions to this problem on non-uniform and uniform meshes with 150 linear finite elements using the standard consistent and lumped mass matrices as well as the averaged mass matrix with reduced dispersion are shown in Figs. 2–5. According to Section 2.2.1 we use the non-uniform mesh with 120 equal finite 1 2 elements of size dx ¼ 120 for the left part of the bar 0 6 x 6 2 (sub2 2 domain p ¼ 1) and 30 equal finite elements of size dx ¼ 30 for the right part of the bar 2 6 x 6 4 (subdomain p ¼ 2). In this case, p co p ¼ const for both subdomains (p ¼ 1; 2). Curves 1 and 2 in dx Fig. 2 correspond to the results obtained in basic computations with the lumped mass matrix and the explicit central-difference method with the time increments equal to the stability limit st (dt ¼ 0:004167 for the non-uniform mesh in (a) and st dt ¼ 0:001666 for the uniform mesh in (b)) and very small time increments dt ¼ 0:0001 for the non-uniform mesh (see curve 2 in Fig. 2(a)). Similar to a homogeneous bar with a uniform mesh, the numerical solution obtained in basic computations by the central-difference method with the time increments equal to the stast bility limit dt ¼ 0:004167 and the lumped mass matrix is close to the exact solution on the non-uniform mesh for the bi-material bar; see curves 1 in Fig. 2. This solution, curves 1, does not have spurious oscillations, does not require the filtering stage and the jumps in velocities in this solution are spread over 3 nodes (for all other points the numerical and exact solutions are the same). If we use a uniform mesh with the lumped mass matrix and the st time increments equal to the stability limit dt ¼ 0:001666 then the numerical solution after basic computations includes large spurious oscillations and is very inaccurate; see curve 2 in Fig. 2(b). The spurious oscillations can be removed by the filtering stage; see curve 3 in Fig. 2(b). Nevertheless, a very fine mesh will be necessary in order to obtain an accurate numerical solution on a uniform mesh. We should also mention that at very small time increments, the numerical solution on a non-uniform mesh includes large spurious oscillations and is also very inaccurate after basic computations; see curve 2 in Fig. 2(a). The spurious oscillations can be removed by the filtering stage (see curve 3 in Fig. 2(a)) and this solution on the non-uniform mesh after the filtering stage is more accurate than that on the uniform mesh; see curves 3 in Fig. 2(a) and (b). Nevertheless, a fine mesh will be necessary in order to obtain an accurate numerical solution on a non-uniform mesh in this case as well. Figs. 3–5 show the numerical results obtained on non-uniform and uniform meshes with 150 linear finite elements using the standard consistent mass matrix as well as the averaged mass matrix with reduced dispersion and the implicit trapezoidal rule with very small time increments dt ¼ 0:0001 in basic computations (in contrast to the lumped mass matrix, smaller time increments yield more accurate results for the considered cases). The numerical solutions after basic computations on the uniform and nonuniform meshes include large spurious oscillations and it is difficult to compare them; see curves 2 and 3 in Figs. 3(a) and 4(a). The results after the filtering stage show that at the same number of degrees of freedom, the non-uniform mesh yields more accurate solutions than the uniform mesh does; see curves 2 and 3 in Figs. 3(b) and 4(b). We can also see that at the same number of degrees of freedom, the averaged mass matrix with reduced

A. Idesman / Composite Structures 117 (2014) 298–308

303

Fig. 2. The velocity distribution along the bar at the observation time T ¼ 98:3. The lumped mass matrix and non-uniform (a) and uniform (b) meshes with 150 linear finite elements are used. Curves 1, 2 and 3 correspond to the results after basic computations and after the filtering stage, respectively. Curves 1 in (a) and (b) are the same and correspond to the non-uniform mesh and the time increments equal to the stability limit for the non-uniform mesh. Curves 2 correspond to the small time increments st st dt ¼ 0:0001 in (a) and the time increments dt ¼ dt ¼ 0:001666 in (b) (where dt is the stability limit for the uniform mesh).

Fig. 3. The velocity distribution along the bar at the observation time T ¼ 98:3. (a) and (b) correspond to the numerical results after basic computations and after the filtering stage. Curve 1 is the exact solution. Curves 2 and 3 correspond to the numerical results obtained on the non-uniform (curve 2) and uniform (curve 3) meshes with 150 linear finite elements. The consistent mass matrix is used.

Fig. 4. The velocity distribution along the bar at the observation time T ¼ 98:3. (a) and (b) correspond to the numerical results after basic computations and after the filtering stage. Curve 1 is the exact solution. Curves 2 and 3 correspond to the numerical results obtained on the non-uniform (curve 2) and uniform (curve 3) meshes with 150 linear finite elements. The averaged mass matrix (with the reduced dispersion) is used.

dispersion yields a much more accurate solution than the consistent matrix does; see curves 2 and 3 in Fig. 5. An example of the convergence of numerical solutions at mesh refinement after the filtering stage is presented in Fig. 6 for the p c averaged mass matrix and non-uniform meshes with dxop ¼ const for both subdomains (p ¼ 1; 2). As can be seen, due to the absence of spurious oscillations, these solutions (curves 2–4 in Fig. 6) converge to the exact solution (curves 1 in Fig. 6). 3.2. Wave propagation in a 1-D bar made of functionally graded materials under impact

Fig. 5. The velocity distribution along the bar at the observation time T ¼ 98:3 after the filtering stage. Curve 1 is the exact solution. Curves 2 and 3 correspond to the numerical results with the averaged mass matrix (with the reduced dispersion) and with the consistent mass matrix. The non-uniform mesh with 150 linear finite elements is used.

In this Section we will consider wave propagation in a 1-D bar of length L ¼ 4 under impact loading (see Fig. 7). The bar is made of a functionally graded material with the following properties: Young’s modulus is chosen to be E ¼ ð1 þ 0:75xÞ2 for 0 6 x 6 4 and the density to be q ¼ 1 for the whole bar where the origin of

304

A. Idesman / Composite Structures 117 (2014) 298–308

Fig. 6. The velocity distribution along the bar at the observation time T ¼ 98:3 after the filtering stage. Curve 1 is the exact solution. Curves 2, 3 and 4 correspond to the numerical results with the averaged mass matrix (with the reduced dispersion) obtained on non-uniform meshes with 150, 300 and 1000 linear finite elements.

Fig. 7. A bar made of a functionally graded material under impact loading.

the Cartesian system is located at the left end of the bar. This corresponds to the following linear variation of the wave velocity: co ¼ ð1 þ 0:75xÞ for 0 6 x 6 4. The following boundary conditions (the same as for the problem in Section 3.1) are applied: the displacement uð0; tÞ ¼ t (which corresponds to the velocity v ð0; tÞ ¼ v 0 ¼ 1) and uð4; tÞ ¼ 0 (which corresponds to the velocity v ð4; tÞ ¼ 0). Initial displacements and velocities are zero; i.e., uðx; 0Þ ¼ v ðx; 0Þ ¼ 0. The observation time is chosen to be T ¼ 99. During this long time the velocity pulse travels within the bar with multiple reflections from each end of the bar. The analytical solution to this problem can be found using the approach developed in [28]. However, this analytical solution is based on the infinite series and may lead to spurious oscillations when a finite number of terms is used in the series; e.g., see [4]. In contrast to this case, the numerical solutions based on the special non-uniform meshes, on the lumped mass matrix and on the explicit central-difference method with the time increments equal to the stability limit yield the numerical results without spurious oscillations and at mesh refinement they converge very fast to the exact solution; see Fig. 8(a). The selection of non-uniform meshes for this problem is described in Section 2.2.2. The numerical solution on the finest mesh with 1000 linear elements will be used as a reference solution; see curve 1 in Fig. 8(a). We should also mention that the filtering stage is not used for the numerical solutions in Fig. 8(a).

In Figs. 9–11 we present the numerical solutions obtained by the two-stage time-integration technique on non-uniform and uniform meshes with 100 linear finite elements using the standard consistent and lumped mass matrices as well as the averaged mass matrix with reduced dispersion. Curves 1 and 2 in Fig. 9 correspond to the results obtained in basic computations with the lumped mass matrix and the explicit central-difference method with the st time increments equal to the stability limit (dt ¼ 0:0046209 for st the non-uniform mesh in (a) and dt ¼ 0:002509 for the uniform mesh in (b)) and very small time increments dt ¼ 0:0001 for the non-uniform mesh (see curve 2 in Fig. 9(a)). Similar to the bi-material bar in Section 3.1, the numerical solution obtained in basic computations by the central-difference method with the time st increments equal to the stability limit dt ¼ 0:0046209 and the lumped mass matrix is close to the exact solution on the nonuniform mesh for the functionally graded bar; see curves 1 in Fig. 9 and curves 1 and 2 in Fig. 8(a) . This solution, curves 1 in Fig. 9, does not have spurious oscillations, does not require the filtering stage and the jumps in velocities in this solution are spread over 3 nodes. If we use a uniform mesh with the lumped mass matrix and the time increments equal to the stability limit st dt ¼ 0:002509 then the numerical solution after basic computations includes large spurious oscillations and is very inaccurate; see curve 2 in Fig. 9(b). The spurious oscillations can be removed by the filtering stage; see curve 3 in Fig. 9(b). Nevertheless, a very fine mesh will be necessary in order to obtain an accurate numerical solution on a uniform mesh. We should also mention that at very small time increments, the numerical solution on the nonuniform mesh includes large spurious oscillations and is also very inaccurate after basic computations; see curve 2 in Fig. 9(a). The spurious oscillations can be removed by the filtering stage (see curve 3 in Fig. 9(a)) and this solution on the non-uniform mesh after the filtering stage is more accurate than that on the uniform mesh; see curves 3 in Fig. 9(a) and (b). Nevertheless, a fine mesh will be necessary in order to obtain an accurate numerical solution on a non-uniform mesh in this case as well. Figs. 10 and 11 show the numerical results obtained on nonuniform and uniform meshes with 100 linear finite elements using the standard consistent mass matrix as well as the averaged mass matrix with reduced dispersion and the implicit trapezoidal rule with very small time increments dt ¼ 0:0001 in basic computations (in contrast to the lumped mass matrix, smaller time increments yield more accurate results for the considered cases). The numerical solutions after basic computations on the uniform and non-uniform meshes include large spurious oscillations and it is difficult to compare them; see curves 2 and 3 in Figs. 10(a) and 11(a). The results after the filtering stage show that at the same number of degrees of freedom, the non-uniform mesh yields much more accurate solutions than the uniform mesh does; see

Fig. 8. The velocity distribution along the bar at the observation time T ¼ 99 after basic computations obtained with the lumped mass matrix and the time increments equal to the stability limit (a), and after the filtering stage obtained with the averaged mass matrix (b); see the text. Curves 1, 2, 3 and 4 in (a) correspond to the numerical results obtained on non-uniform meshes with 1000, 100, 50 and 20 linear finite elements. Curve 1 in (b) is the reference solution and is the same as curve 1 in (a). Curves 2, 3 and 4 in (b) correspond to the numerical results obtained on non-uniform meshes with 100, 300 and 1000 linear finite elements.

305

A. Idesman / Composite Structures 117 (2014) 298–308

Fig. 9. The velocity distribution along the bar at the observation time T ¼ 99. The lumped mass matrix and non-uniform (a) and uniform (b) meshes with 100 linear finite elements are used. Curves 1, 2 and 3 correspond to the results after basic computations and after the filtering stage, respectively. Curves 1 in (a) and (b) are the same and correspond to the non-uniform mesh and the time increments equal to the stability limit for the non-uniform mesh. Curves 2 correspond to the small time increments st st dt ¼ 0:0001 in (a) and the time increments dt ¼ dt ¼ 0:002509 in (b) (where dt is the stability limit for the uniform mesh).

Fig. 10. The velocity distribution along the bar at the observation time T ¼ 99. (a) and (b) correspond to the numerical results after basic computations and after the filtering stage. Curve 1 is the reference solution. Curves 2 and 3 correspond to the numerical results obtained on the non-uniform (curve 2) and uniform (curve 3) meshes with 100 linear finite elements. The consistent mass matrix is used.

Fig. 11. The velocity distribution along the bar at the observation time T ¼ 99. (a) and (b) correspond to the numerical results after basic computations and after the filtering stage. Curve 1 is the reference solution. Curves 2 and 3 correspond to the numerical results obtained on the non-uniform (curve 2) and uniform (curve 3) meshes with 100 linear finite elements. The averaged mass matrix (with the reduced dispersion) is used.

curves 2 and 3 in Fig. 10(b) for the consistent mass matrix and in Fig. 11(b) for the averaged mass matrix. We should mention that the two-stage time-integration technique yields non-oscillatory and convergent numerical solutions for any type of the mass matrix. However, except the case of non-uniform meshes, the lumped mass matrix and the explicit central difference method with the time increments equal to the stability limit, in all other cases the averaged mass matrix yields much more accurate results at the same number of degrees of freedom, see Figs. 9–11. The convergence of numerical solutions at mesh refinement after the filtering stage is presented in Fig. 8(b) for the averaged mass matrix and non-uniform meshes. As can be seen, due to the absence of spurious oscillations after basic computations, these solutions (curves 2–4 in Fig. 8(b)) converge to the reference solution (curve 1 in Fig. 8(b)). We should also mention that without the filtering stage the numerical results (even for the finite elements with reduced dispersion) diverge at

mesh refinement due to large spurious high-frequency oscillations, see Fig. 12. 3.3. Wave propagation in a bi-material cylinder under impact Here we consider an axisymmetric problem of wave propagation in a bi-material cylinder of length L ¼ 4 and radius R ¼ 1 under impact loading; see Fig. 13(a). The cylinder is made of two materials with the following properties: Young’s modulus is chosen to be E ¼ 1 for 0 6 z 6 2 and E ¼ 16 for 2 6 z 6 4; Poisson’s ration to be m ¼ 0:3 and the density to be q ¼ 1 for the whole cylinder where the origin of the Cartesian system is located at the left end of the cylinder at point A and AD is the axis of revolution. These material properties correspond to the following ratios of the wave velocities in two subdomains:

csubdomain2 m csubdomain1 m

¼ 4 where m ¼ 1 and m ¼ 2

correspond to the compressional and shear waves, respectively; subdomains 1 and 2 correspond to 0 6 z 6 2 and 2 6 z 6 4; see

306

A. Idesman / Composite Structures 117 (2014) 298–308

Fig. 12. The velocity distribution along the bar after basic computations at the observation time T ¼ 99. Curve 1 is the reference solution. Curves 2 in (a), (b) and (c) correspond to the numerical results obtained on non-uniform meshes with 100, 300 and 1000 linear finite elements, respectively. The averaged mass matrix (with the reduced dispersion) is used.

(a)

(b) Fig. 13. A bi-material cylinder under impact loading (a). AD is the axis of revolution. An example of a non-uniform mesh (40 and 10 uniform square elements along the radius are used in regions I and III, respectively).

Fig. 13(a). The following boundary conditions are applied: the axial displacement uz ðr; z ¼ 0; tÞ ¼ t along AB (which corresponds to the axial velocity v z ðr; z ¼ 0; tÞ ¼ v 0 ¼ 1) and zero tractive forces along the tangential direction of AB; zero tractive forces along BC and CD. Initial displacements and velocities are zero; i.e., ur ðr; z; 0Þ ¼ uz ðr; z; 0Þ ¼ v r ðr; z; 0Þ ¼ v z ðr; z; 0Þ ¼ 0. The observation time is chosen to be T ¼ 6. During this time the velocity pulse travels within the cylinder with several reflections from each end of the cylinder and from the interface between two materials at z ¼ 2. As shown in Section 3.1, the linear elements with reduced dispersion are much more accurate than the standard linear elements in the 1-D case. The same is valid for the axisymmetric problem

considered here. Therefore, for the axisymmetric problem we will present the numerical results obtained on non-uniform and uniform meshes with 4-node bilinear finite elements and the averaged mass matrix with reduced dispersion. In contrast to the 1-D problem in Section 3.1, in the 2-D case we cannot generate non-uniform meshes consisting of different uniform square elements for subdomains 1 and 2 (in order to meet the condition that the element size is proportional to the wave velocity). Therefore, we will use the following approach for the generation of non-uniform meshes: we will use uniform meshes with square elements for subdomain 1 and for a large part of subdomain 2 where the size of the element in the latter is four times as large

A. Idesman / Composite Structures 117 (2014) 298–308

as that in subdomain 1. Then, for the transition from a fine mesh in subdomain 1 to a coarse mesh in subdomain 2 we use a small transition region with a non-uniform mesh. An example of such a mesh is shown in Fig. 13(b) with 40 and 10 square elements along the radius in regions I and III. We should also mention that the transition region II is located in subdomain 2. Therefore, the use of small elements in the transition region (smaller than those in region III) does not lead to the degradation of accuracy. The two-stage time integration technique with small time increments in basic computations is also applied. Fig. 14 shows the distribution of the axial velocity in the cylinder at the observation time T ¼ 6. The results in Figs. 14(a) and (b) correspond to non-uniform meshes with 160 (a) (320 (b)) and 40 (a) (80 (b)) square finite elements along the radius in AB and CD, respectively (similar to the mesh shown in Fig. 13(b)). As can be seen from Fig. 14(a) and (b) the numerical results without the filtering stage

307

include large spurious oscillations (see curves 2) which can be efficiently removed by the filtering stage (see curves 1). We should emphasize that the results without the filtering stage diverge at mesh refinement due to the increase in the amplitudes of spurious highfrequency oscillations for finer meshes; see curves 2 in Fig. 14(a) and (b). We have also solved the same problem using uniform meshes with 160 and 320 linear elements along the radius for the entire cylinder and have obtained approximately the same results after the filtering stage as curves 1 in Fig. 14(a) and (b). However, the use of the non-uniform meshes allows the reduction in the number of degrees of freedom by a factor of 1.84 at the same accuracy. Fig. 14(c) and (d) shows the convergence of the axial velocity distribution along the axis of revolution AD and along the external surface BC at mesh refinement where curves 1, 2, 3 and 4 correspond to non-uniform meshes with 40; 80; 160 and 320 square elements along BC (and 10; 20; 40 and 80 square elements along CD), respectively. As can be seen from Fig. 14(c) and (d), the variation of the axial velocity is more complicated for subdomain I. We should also mention that despite the uniform impact loading along the radius at the left end of the cylinder, the distribution of the axial velocity along the radius is very non-uniform for both subdomains; compare the distribution of the axial velocity along AD and BC in Fig. 14(c) and (d). 4. Concluding remarks One of the main issues in the finite element modeling of wave propagation problems (especially at impact loading) is the appearance of large spurious high-frequency oscillations that may destroy numerical solutions. In our papers [14–20] we have developed a new finite element approach for the accurate modeling of wave propagation in homogeneous materials. In this paper we have extended this technique to inhomogeneous materials including composite and functionally graded materials and we have shown the main differences for the modeling of wave propagation in inhomogeneous and homogeneous materials. Below we will summarize the main findings of the paper.

Fig. 14. The distribution of the axial velocity V z along the axis of revolution AD (a), (b) and (c) along the external surface BC (d) at the observation time T ¼ 6. Curves 1 and 2 in (a) and (b) correspond to the results on non-uniform meshes with 56222 nodes (a) and 220594 nodes (b) with (curves 1) and without (curves 2) the filtering stage. Curves 1, 2, 3 and 4 in (c) and (d) correspond to the results with the filtering stage and non-uniform meshes with 3633 nodes, 14261 nodes, 56222 nodes and 220594 nodes, respectively.

1. For the special design of finite element meshes (see Section 2.2), we can obtain very accurate results for wave propagation problems in inhomogeneous materials if we use the standard linear elements with the lumped mass matrix and the explicit central difference method with the time increments equal to the stability limit. In this case, there are no spurious oscillations in the numerical solutions after basic computations (the filtering stage is not required) and a discontinuity in a solution (if presented) is spread over three nodes. In contrast to the finite element formulations with graded finite elements (e.g., see [4,6,10,12]), we use piecewise constant variation of materials properties (they are constant within any finite element). Nevertheless, even for functionally graded materials the numerical results converge very fast to the exact solution at mesh refinement. It is interesting to mention that in the 1-D case the analytical solution to wave propagation in functionally graded materials can be obtain by the approach developed in [28]. However, this analytical solution is based on the infinite series and may lead to spurious oscillations due to Gibbs phenomena when a finite number of terms is used in the series; e.g., see [4]. In contrast to this case, the numerical solutions based on the special non-uniform meshes with standard linear finite elements, on the lumped mass matrix and on the explicit central-difference method with the time increments equal to the stability limit do not include spurious oscillations at any number of degrees of freedom. They can be used as reference solutions for testing computer codes as well as for the analysis of wave propagation phenomena under highfrequency and impact loadings.

308

A. Idesman / Composite Structures 117 (2014) 298–308

2. Except one case described above, in all other cases (e.g., uniform meshes or the non-lumped mass matrix or smaller time increments or other time-integration methods or high-order finite elements or the 2-D and 3-D problems etc.) numerical solutions obtained by existing methods may include large spurious high-frequency oscillations especially under impact loading. However, the new numerical approach based on the two-stage time-integration technique (with basic computations and the filtering stage) quantifies and removes the spurious oscillations for these cases. For the first time we have obtained accurate finite element solutions of wave propagation problems in composite and functionally graded materials without spurious oscillations and without the interaction between user and computer code. 3. In many papers on the simulation of the dynamic response of structures made of inhomogeneous materials (e.g., see [1–12] for composite and functionally graded materials) uniform meshes or meshes independent of material properties are used with the standard approaches with explicit or implicit timeintegration methods. In the paper we have shown that despite uniform meshes are justified for homogeneous materials, nonuniform meshes with the constant ratios of the local wave velocity to the size of the finite element yield much more accurate results than uniform meshes do. 4. We have shown that for inhomogeneous materials (similar to homogeneous materials), the finite element formulation with reduced dispersion (based on the averaged mass matrix and the linear elements) yields much more accurate results than the standard formulations with the lumped or consistent mass matrices do. There is one exception for the case of the lumped mass matrix with the non-uniform mesh with linear elements and the explicit central difference method with the time increments equal to the stability limit. However, this exception is valid in the 1-D case only. In our paper [19] we have shown that for homogeneous materials in the multidimensional case, the formulation with reduced dispersion is much more efficient than other formulations. We should also mention that in this paper we use the finite element formulation with reduced dispersion along with the implicit time-integration method. However, similar results can be obtained for the finite element formulation with reduced dispersion used with explicit timeintegration methods; see our paper [19]. 5. In our previous papers [15,16,19,20] we have shown that the size of time increments at the filtering stage of the two-stage time-integration technique (see Eqs. (7)–(9)) defines the range of actual frequencies used in numerical solutions and can serve as a quantitative measure for the comparison and the prediction of the accuracy and the computational efficiency of different space-discretization techniques used with different meshes. The smaller size of time increments at the filtering stage calculated according to Eqs. (7)–(9) corresponds to more accurate numerical solutions; e.g., see our papers [15,16,19,20]. 6. As shown in our paper [14], the new approach can be equally applied to all elastodynamics problems including structural dynamics and wave propagation problems (independent of the initial and boundary conditions).

Acknowledgments The research has been supported in part by the Air Force Office of Scientific Research (contract FA9550-12-1-0324) and by Texas Tech University.

References [1] Aksoy H, Senocak E. Wave propagation in functionally graded and layered materials. Finite Elem Anal Des 2009;45(12):876–91. [2] Banks-Sills L, Eliasi R, Berlin Y. Modeling of functionally graded materials in dynamic analyses. Compos Part B : Eng 2002;33(1):7–15. [3] Chakraborty A, Gopalakrishnan S. Wave propagation in inhomogeneous layered media: solution of forward and inverse problems. Acta Mech 2004;169(1-4):153–85. [4] Santare MH, Thamburaj P, Gazonas GA. The use of graded finite elements in the study of elastic wave propagation in continuously nonhomogeneous materials. Int J Solids Struct 2003;40(21):5621–34. [5] Jung W-Y, Han S-C. Transient analysis of fgm and laminated composite structures using a refined 8-node ans shell element. Compos Part B : Eng 2014; 56:372–83. [6] Kim J-H, Paulino G. Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials. J Appl Mech, Trans ASME 2002;69(4): 502–14. [7] Mahapatra DR, Singhal A, Gopalakrishnan S. A higher-order finite waveguide model for spectral analysis of composite structures. Comput Methods Appl Mech Eng 2006;195(9-12):1116–35. [8] Marur P, Tippur H. Dynamic response of bimaterial and graded interface cracks under impact loading. Int J Fract 2000;103(1):95–109. [9] Oyekoya O, Mba D, El-Zafrany A. Buckling and vibration analysis of functionally graded composite structures using the finite element method. Compos Struct 2009;89(1):134–42. [10] Park K, Paulino GH. Parallel computing of wave propagation in three dimensional functionally graded media. Mech Res Commun 2011;38(6): 431–6. [11] Yas M, Shakeri M, Heshmati M, Mohammadi S. Layer-wise finite element analysis of functionally graded cylindrical shell under dynamic load. J Mech Sci Technol 2011;25(3):597–604. [12] Zhang ZJ, Paulino GH. Wave propagation and dynamic analysis of smoothly graded heterogeneous continua using graded finite elements. Int J Solids Struct 2007;44(11-12):3601–26. [13] Hui D, Dutta PK. A new concept of shock mitigation by impedance-graded materials. Compos Part B : Eng 2011;42(8):2181–4. [14] Idesman AV. Accurate time integration of linear elastodynamics problems. Comput Model Eng Sci 2011;71(2):111–48. [15] Idesman A, Schmidt M, Foley JR. Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error. Comput Mech 2011;47:555–72. [16] Idesman A, Schmidt M, Foley JR. Accurate 3-d finite element simulation of elastic wave propagation with the combination of explicit and implicit timeintegration methods. Wave Motion 2011;48:625–33. [17] Idesman AV, Samajder H, Aulisa E, Seshaiyer P. Benchmark problems for wave propagation in elastic materials. Comput Mech 2009;43(6):797–814. [18] Idesman AV, Subramanian K, Schmidt M, Foley JR, Tu Y, Sierakowski RL. Finite element simulation of wave propagation in an axisymmetric bar. J Sound Vib 2010;329:2851–72. [19] Idesman A, Pham D. Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error. Comput Methods Appl Mech Eng 2014;271:86–108. [20] Idesman A, Pham D, Foley JR, Schmidt M. Accurate solutions of wave propagation problems under impact loading by the standard, spectral and isogeometric high-order finite elements. comparative study of accuracy of different space-discretization techniques. Finite Elem Anal Des 2014;88:67–89. [21] Idesman A, Mates SP. Accurate finite element simulation and experimental study of elastic wave propagation in a long cylinder under impact loading. Int J Impact Eng 2014;71:1–16. [22] Cherukuri HP. Dispersion analysis of numerical approximations to plane wave motions in an isotropic elastic solid. Comput Mech 2000;25(4):317–28. [23] Dauksher W, Emery AF. Solution of elastostatic and elastodynamic problems with chebyshev spectral finite elements. Comput Methods Appl Mech Eng 2000;188(1):217–33. [24] Guddati MN, Yue B. Modified integration rules for reducing dispersion error in finite element method. Comput Methods Appl Mech Eng 2004;193:275–87. [25] Krenk S. Dispersion-corrected explicit integration of the wave equation. Comput Methods Appl Mech Eng 2001;191:975–87. [26] Yue B, Guddati MN. Dispersion-reducing finite elements for transient acoustics. J Acoust Soc Am 2005;118(4):2132–41. [27] Hughes T, Reali A, Sangalli G. Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of p-method finite elements with k-method nurbs. Comput Methods Appl Mech Eng 2008;197(49-50):4104–24. [28] Chiu T-C, Erdogan F. One-dimensional wave propagation in a functionally graded elastic medium. J Sound Vib 1999;222(3):453–87.