A mixed Ritz-DQ method for forced vibration of functionally graded beams carrying moving loads

Composite Structures 92 (2010) 2497–2511

Contents lists available at ScienceDirect

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

A mixed Ritz-DQ method for forced vibration of functionally graded beams carrying moving loads S.M.R. Khalili a,b,*, A.A. Jafari a, S.A. Eftekhari a a b

Centre of Excellence for Research in Advanced Materials and Structures, Faculty of Mechanical Engineering, K.N. Toosi University of Technology, Tehran, Iran Faculty of Engineering, Kingston University, London, UK

a r t i c l e

i n f o

Article history: Available online 26 February 2010 Keywords: DQM Rayleigh–Ritz method FG beam Moving load Inertia effect Forced vibration

a b s t r a c t A mixed method is presented to study the dynamic behavior of functionally graded (FG) beams subjected to moving loads. The theoretical formulations are based on Euler–Bernoulli beam theory, and the governing equations of motion of the system are derived using the Lagrange equations. The Rayleigh–Ritz method is employed to discretize the spatial partial derivatives and a step-by-step differential quadrature method (DQM) is used for the discretization of temporal derivatives. It is shown that the proposed mixed method is very efficient and reliable. Also, compared to the single-step methods such as the Newmark and Wilson methods, the DQM gives better accuracy using larger time step sizes for the cases considered. Moreover, effects of material properties of the FG beam and inertia of the moving load on the dynamic behavior of the system are investigated and analyzed. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction The study of behavior of functionally graded materials (FGMs) has been an interesting topic of considerable focus during the past decades. The intensity and rapid growth of research into this class of materials is actually due to their continuously varying material properties, which gives great advantages over the conventional homogeneous and layered materials. For example, FGMs consisting of metallic and ceramic components used in dynamic systems can improve both the mechanical and thermo-mechanical performances of the system. On the other hand, the cracking and delamination phenomenon which are often observed in conventional multi-layered systems are avoided due to the smooth transition between the properties of the components in FGMs [1–3]. Therefore, it is of great importance to analyze the behavior of FG structures. Due to the above-mentioned favorable features, many studies have been conducted on the static and dynamic behavior of FG structures. The elasticity solutions of transversely and thermally loaded FG beams and also sandwich beams with FG cores were obtained by Sankar and his co-workers [4–6]. In their studies, the thermo mechanical properties of FGM were all assumed to vary exponentially through the thickness of the beam. Zhu and Sankar [7] solved the two-dimensional elasticity equations for the FG beam

* Corresponding author. Address: Centre of Excellence for Research in Advanced Materials and Structures, Faculty of Mechanical Engineering, K.N. Toosi University of Technology, Pardis St., Molasadra Ave., Vanak Sq., Tehran, Iran. Tel.: +98 2184063208; fax: +98 2188674748. E-mail address: [email protected] (S.M.R. Khalili). 0263-8223/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2010.02.012

subjected to transverse loads using a mixed Fourier series – Galerkin method. Shi and Chen [8] studied the problem of a functionally graded piezoelectric cantilever beam subjected to different loadings. A pair of stress and induction functions was proposed in the form of polynomials and a set of analytical solutions for the beam subjected to different loadings was obtained. Nirmala et al. [9] derived analytical solutions for the thermo-elastic stresses in a threelayer composite beam with a middle FG layer. Using the meshless local Petrov–Galerkin method, Ching and Yen [10] presented numerical solutions for two-dimensional FG solids such as circular cylinders and simply supported beams subjected to either mechanical or thermal loads. Later, they used the meshless technique to obtain the transient thermo-elastic responses of FG beams under a non-uniform connective heat supply [11]. The free vibration of orthotropic FG beams with various end conditions was studied by Lu and Chen [12]. Wu et al. [13] obtained the closed-form natural frequencies of a simply supported FG beam with mass density and Young’s modulus being polynomial functions of the axial coordinate. Aydogdu and Taskin [14] and Yang and Chen [15] analyzed the free vibration of FG beams. Kapuria et al. [16] presented a finite element model for static and free vibration responses of layered FG beams using a third-order theory and its experimental validation. A plane elasticity solution for a cantilever FG beam subjected to different static loads by means of the semi-inverse method was derived by Zhong and Yu [17]. Two-dimensional elasticity solutions for bending and free vibration of FG beams resting on Winkler–Pasternak elastic foundation were obtained by Ying et al. [18]. Li [19] proposed a new unified approach to investigate the static and free vibration behavior of Euler–Bernoulli and Timoshenko beams.

2498

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

Nomenclature m nT Aij Bij Et Eb Eratio

qt qb qratio

number of sampling time points per DQM time element number of DQM time elements elements of first-order DQM weighting coefficient matrix elements of second-order DQM weighting coefficient matrix modulus of elasticity (Young’s modulus) at the top surface of the FG beam modulus of elasticity (Young’s modulus) at the bottom surface of the FG beam ratio of Young’s modulus at the top and bottom surfaces of the FG beam (=Et/Eb) density at the top surface of the beam density at the bottom surface of the beam ratio of density at the top and bottom surfaces of the FG beam (=qt/qb)

Based on the high-order sandwich panel theory, Rahmani et al. [20] developed a new model for free vibration analysis of sandwich beams with syntactic foam as a FG flexible core. From the review of the literature, it is found that most of the researchers are interested in the free vibration analysis of FG structures, but few of them have paid attention to the forced vibration analysis of FG structures. On the other hand, to the authors’ best knowledge, the information regarding the forced vibration responses of FG structures due to the moving loads is rare and this is the reason why this paper aims its study. In the literature, one may refer to the work done by Yang et al. [21] or the recent study by Sßimsßek and Kocatürk [22]. In Ref. [21], the free and forced vibrations of cracked FG beams subjected to an axial force and a moving load were studied by using the modal technique. In the recent paper [22], the free and forced vibrations of a simply supported FG beam subjected to a concentrated moving harmonic load were analyzed. The transverse and axial displacements of the simply supported beam were first approximated by series of polynomial functions and by using the Lagrange equations the governing equations of motion of the FG beam were obtained in matrix form. The Newmark time integration scheme was then employed to solve the resulting dynamic equations. In both the papers [21,22], the effect of inertia of the moving load was ignored and the moving load was modeled as a simple moving force model (a force with constant [21] or harmonic magnitude [22] which moved along the FG beam). Besides, in Ref. [21] the material properties of the beam were only assumed to vary exponentially in the beam thickness direction and also the results of Ref. [22] were only given for simply supported FG beams. The differential quadrature method (DQM), which was first introduced by Bellman et al. [23,24] in the early 1970s, is an alternative numerical solution technique for initial and/or boundary problems in engineering and mathematics. Its central idea is to approximate the derivative of a function at a discrete point with respect to a coordinate direction using a linear sum of all the function values at all discrete points along that direction [25,26]. Compared to the low-order methods, such as the finite element and finite difference methods, the DQ method can generate numerical results with high-order of accuracy by using a considerably smaller number of discrete points and therefore requiring relatively little computational effort. Another particular advantage of the DQ method lies in its ease of use and implementation. The DQ method was initially proposed for the solution of initialvalue problems [23], but little attention has been paid on this issue until recent years. Using the Hermite interpolation functions as the trial functions, Wu and Liu [27,28] proposed a generalized differen-

x horizontal coordinate z vertical coordinate t time k power-law exponent u axial displacement of the beam w vertical displacement of the beam /(x), w(x) Rayleigh–Ritz trial/test functions M mass of the moving load v velocity of the moving point load g gravitational acceleration [m*], [c*], [k*] instantaneous mass, damping, and stiffness matrices due to the inertia of the moving load [M], [C], [K] structural mass, damping, and stiffness matrices

tial quadrature rule to solve linear and nonlinear initial-value problems. Shu and Yao [29] proposed a block Marching technique with DQ discretization for the solution of time-dependent problems. The time span was first divided into several blocks and the quadrature rule was then applied in each of the blocks. This technique reduces considerably the computational cost, since a smaller system of algebraic equations should be solved at each block (step). Tanaka and Chen [30,31] proposed a combined application of dual reciprocity boundary element method (DRBEM) and DQ method to transient diffusion and elastodynamic problems. In a series of papers on the solution of initial-value problems using the DQ method, Fung [32–34] introduced a modified DQ method to incorporate initial conditions. He also discussed in details, the stability property of the DQ method. It was emphasized that the accuracy and stability of the DQ method are dictated by the choice of sampling time points and also by the ways in which the initial conditions are incorporated into the solution process. He also showed that unconditionally stable algorithms using the DQ method can be obtained, if the initial conditions are incorporated by modifying the weighting coefficient matrices and if the number of sampling time points is small. Based on the construction of cubic cardinal spline functions, Zhong and Lan [35] developed a DQ time integration scheme and applied it to the solution of dynamic systems governed by Duffing-type nonlinear differential equations. Eftekhari et al. [36–40] proposed a coupled finite element – differential quadrature method for the solution of moving load class of problem. It has been claimed that the proposed mixed method is very efficient for the solution of time-dependent problems. Liu and Wang [41] presented an assessment of the DQ time integration scheme for nonlinear dynamic equations. It was shown that the DQ time integration scheme is reliable, computationally efficient and also suitable for time integrations over long time duration. But care should be taken in choosing a time step when applying the DQ method to nonlinear systems. In a series of papers, Civalek proposed a number of mixed methodologies for the nonlinear dynamic analysis of rectangular plates and multi-degree-of-freedom (MDOF) systems [42–44]. In this study a new mixed method, which combines the Rayleigh–Ritz method and the differential quadrature method (DQM), is proposed for forced vibration analysis of FG beams subjected to moving loads. The governing equations of motion of the FG beam are derived by using the Lagrange equations under the assumption of the Euler–Bernoulli beam theory. The Rayleigh–Ritz method is first employed to discretize the spatial partial derivatives, and the DQM is then applied to analogize the temporal derivatives. The resulting system of algebraic equations can be easily

2499

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

solved using various direct (such as the Gaussian elimination method) or iterative (such as the Gauss–Seidel method) methods. To obtain the element of structural mass and stiffness matrices, it is assumed that the material properties of the FG beam vary continuously in the thickness direction according to the exponential law and power-law functions. The moving load is modeled as a moving mass load to account the effect of inertia of the moving load. To the best of the authors’ knowledge, this is the first attempt on the study of forced vibration of FG beams subjected to concentrated moving masses. Since the velocity of moving load, the inertia of the moving load and the ratio of material properties of the FG beam are the key factors affecting the dynamic behavior of the beam, the influence of these parameters are investigated. The dynamic equations of motion of the beam are also solved by using implicit time step methods such as the Newmark and Wilson methods. It is shown that the DQ time integration scheme can produce much better accuracy than the Newmark and Wilson methods using much larger time step sizes and therefore needing relatively little computational cost.

The physical system analyzed here, is a functionally graded Euler–Bernoulli beam of length L, width b and thickness h, supporting a mass moving at a velocity v, as shown in Fig. 1. The mass enters the beam at t = 0 and leaves the beam at t = T. The velocity of the moving mass is assumed to be constant. In this study, the material properties (i.e., Young’s modulus E and density q) of the FG beam are assumed to vary through the thickness according to exponential and power-low functions. The exponential law, which is more common in fracture studies of FGM [45,46] is given by:

ð1Þ

The power-law, which was first introduced by Wakashima et al. [47] is given by:

 k 1 z PðzÞ ¼ ðP t  P b Þ þ Pb þ 2 h

S¼ K¼

1 2

Z

1 2

ð6Þ

qðzÞðu_ 2 þ w_ 2 ÞdA dx

ð7Þ

Z

L

A

0

Substituting Eqs. (3)–(5) in Eqs. (6) and (7) leads to:

S¼ K¼

1 2

Z

1 2

L

0

Z

  Axx u20;x  2Bxx u0;x w0;xx þ Dxx w20;xx dx

ð8Þ

  IA u20;t  2IB u0;t w0;xt þ IA w20;t þ ID w20;xt dx

ð9Þ

L

0

where

Z ðAxx ; Bxx ; Dxx Þ ¼ EðzÞð1; z; z2 Þdz A Z qðzÞð1; z; z2 Þdz ðIA ; IB ; ID Þ ¼

ð10Þ ð11Þ

A

W¼

Z

ð3Þ

L

0

€ 0 ðx; tÞjx¼v t Þdðx  v tÞw0 ðx; tÞdx ðMg þ M w

€ 0 ðv t; tÞÞw0 ðv t; tÞ W ¼ Mðg þ w

ð13Þ

The n-parameter Rayleigh–Ritz solutions for the problem are of the following forms [48]:

w0 ðx; tÞ ¼

n X

wj ðtÞ/j ðxÞ

ð14Þ

uj ðtÞwj ðxÞ

ð15Þ

j¼1 n X j¼1

where wj and uj (j = 1, 2, . . ., n) are Ritz coefficients and /j(x) and wj(x) are Ritz approximation/test functions. The Rayleigh–Ritz test functions should satisfy at least geometric boundary conditions (see Appendix A for more details). If the boundary conditions are non-homogeneous, then some terms should be added to Eqs. (14) and (15) to satisfy these conditions. Now, from Eq. (14) and successful differentiation of this equation with respect to time, it is obtained that:

where u0 and w0 are the axial and transverse displacements of the beam in the reference plane (the coordinates (x, z) shown in Fig. 1), respectively, and (),x represents the partial differentiations with respect to x. Using Eq. (3), the linear strain and the resultant stress can be written as:

w0 ðv t; tÞ ¼

exx ¼ u;x ¼ u0;x  zw0;xx rxx ¼ EðzÞexx

€ 0 ðv t; tÞ ¼ w

ð4Þ

ð12Þ

where the over-double-dot denotes double differentiation with respect to time, M is the mass of moving load, v is the velocity of moving load, g is the gravitational acceleration (9.81 m/s2), and d() denotes the Dirac delta-function. The first and second terms of Eq. (12) are the potential energy due to the weight and inertia of the moving load, respectively. From Eq. (12), it is obtained that:

ð2Þ

P(z) denotes a typical material property (E, q). Pt and Pb denote the values of the variables at the top and bottom of the FG beam, respectively, and k is a parameter describing the material variation profile through the thickness of the FG beam. Based on the Euler–Bernoulli beam theory, the axial and transverse displacements of the beam are expressed as

n X

wj ðtÞ/j ðv tÞ

ð16Þ

j¼1

_ 0 ðv t; tÞ ¼ w

n X _ j ðtÞ/j ðv tÞ þ v wj ðtÞ/j;x ðv tÞÞ ðw

ð17Þ

j¼1 n X € j ðtÞ/j ðv tÞ þ 2v w _ j ðtÞ/j;x ðv tÞ þ v 2 wj ðtÞ/j;xx ðv tÞÞ ðw j¼1

ð5Þ

z

rxx exx dA dx

A

0

Z

u0 ðx; tÞ ¼

uðx; z; tÞ ¼ u0 ðx; tÞ  zw0;x ; wðx; z; tÞ ¼ w0 ðx; tÞ

Z

L

The potential energy of the moving load can be expressed as:

2. Theory and Rayleigh–Ritz formulation

  1 Pt PðzÞ ¼ P t expðd0 ð1  2z=hÞÞ; d0 ¼ ln 2 Pb

The strain energy (S) and the kinetic energy (K) of the beam are given by:

ð18Þ Note that the velocity of the moving mass (i.e., v) is constant and

v M x

/j;x ðv tÞ ¼

 d/j ðxÞ ; dx x¼v t

/j;xx ðv tÞ ¼

 2 d /j ðxÞ  2 dx 

ð19Þ

x¼v t

The Lagrangian of the present problem is: Fig. 1. Simply supported FG beam with a moving mass.

I ¼ K  ðS þ WÞ

ð20Þ

2500

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

Using the Lagrange equation, the governing equation of motion of the system can be obtained as:

"

#



  € _ fwg fwg ½cðtÞ ½0 þ €g _ fu fug ½0 ½0 ½M 22  ½M 21  " #     11 12 fwg ff ðtÞg ½K  þ ½kðtÞ ½K  ¼ þ fug f0g ½K 22  ½K 21 

½M 11  þ ½mðtÞ

½M 12 

where [0] is the n  n zero matrix, fwg ¼ ½ w1 fug ¼ ½ u1 u2 . . . un T and

M 11 ij

¼ IA

Z

L

/i /j dx þ ID

Z

0

M 22 ij ¼ IA

wn T ;

/i;x /j;x dx

ð22Þ ð23Þ

€f ðt Þ ¼ i

wi /j;x dx

ð24Þ

wi wj dx

0

ð25Þ

L

/i;xx /j;xx dx

ð26Þ

0

Z

L

/i;xx wj;x dx

ð27Þ

wi;x /j;xx dx

ð28Þ

0

Z

L

0 L

Z 0

m X

Aij f ðt j Þ

ð37Þ

Bij f ðtj Þ

ð38Þ

j¼1 m X j¼1

L

Z

K 21 ij ¼ Bxx

Let f(t) be a solution of a differential equation and t1 = 0, t2, t3, . . ., tm = T be a set of sample points in the time domain. According to DQM, the first- and second-order derivatives of the function f(t) at any sample time points can be approximated by the following formulation [23]:

f_ ðt i Þ ¼

0 L

K 12 ij ¼ Bxx

K 22 ij ¼ Axx

...

L

/i;x wj dx

0

Z

K 11 ij ¼ Dxx

ð21Þ

L

Z

M 21 ij ¼ I B

3. Differential quadrature method (DQM)

0

Z

M 12 ij ¼ I B

w2

are neglected, the problem is referred to moving force problem and otherwise is referred to moving mass problem.

wi;x wj;x dx

ð29Þ

fi ðtÞ ¼ Mg/i ðv tÞ

ð30Þ

The time-dependent instantaneous matrices [m(t)]*, [c(t)]*, and [k(t)]* (due to the inertia of the moving load) are given by:

mij ðtÞ ¼ M/i ðv tÞ/j ðv tÞ

ð31Þ

cij ðtÞ ¼ 2M v /i ðv tÞ/j;x ðv tÞ

ð32Þ

 kij ðtÞ

ð33Þ

2

¼ M v /i ðv tÞ/j;xx ðv tÞ

Eq. (21) can also be written in the compact form:

€ðtÞg þ ½CðtÞfyðtÞg _ ½MðtÞfy þ ½KðtÞfyðtÞg ¼ fFðtÞg

ð34Þ

where m is the number of sample points in the time domain, f(tj) represents the functional value at a sample time point tj, f_ ðt i Þ and €f ðt Þ indicate the first- and second-order derivatives of f(t) at a time i point ti, Aij and Bij are the weighting coefficients of the first- and second-order derivatives, respectively. The weighting coefficients can be determined by the functional approximations in the time domain. Using the Lagrange interpolation polynomials as the approximating functions (say test functions), Quan and Chang [49] obtained the following algebraic formulations to compute the first-order weighting coefficients:

8 Mð1Þ ðt i Þ > > < ðti tk ÞMð1Þ ðtk Þ i – k; i; k ¼ 1; 2; . . . ; m Aik ¼ m P > > Aij i ¼ k; i ¼ 1; 2; . . . ; m :

ð39Þ

j¼1;j–i

(1)

where M (t) is defined as:

Mð1Þ ðt i Þ ¼

m Y

ðti  tj Þ

ð40Þ

j¼1;j – i

The weighting coefficients of the second-order derivative can be obtained from the following recurrence relationship [25]:

8 h i 1 > i – k; i; k ¼ 1; 2; . . . ; m > < 2Aik Aii  ti tk Bik ¼ m P > > Bij i ¼ k; i ¼ 1; 2; . . . ; m :

ð41Þ

j¼1;j–i

where the transient mass, damping, and stiffness matrices are of order 2n  2n and the displacement and force vectors are of order 2n  1. Mathematically, Eq. (34) is a system of coupled ordinary differential equations of second-order in time with time-dependent coefficient matrices, which can be solved using various explicit or implicit time step methods. In this study, Eq. (34) is solved by using the DQ time integration scheme. The details of this scheme is described in the next section. Note that when the traveling mass leaves the beam, it can be obtained that:

½m ¼ ½c ¼ ½k ¼ ½0;

fFg ¼ f0g

ð35Þ

Therefore, the governing equations of motion of the system read to be:

€ðtÞg þ ½KfyðtÞg ¼ f0g ½Mfy

ð36Þ

Also, note that when the effects of inertia of the moving load are neglected, the mass and stiffness matrices are no longer timedependent. In this paper, when effects of inertia of the moving load

One of the key factors in the accuracy and rate of convergence of the DQ solutions is the choice of sampling points. It is well known that the equally spaced sampling points are not very desirable [26]. It has been suggested that non-uniformly spaced sampling points can generate more accurate solutions. The zeros of some of the orthogonal polynomials are commonly adopted as the sampling points. In this work, the Chebyshev–Gauss–Lobatto sample points are used for calculation of weighting coefficients. These points are given by [26]:

 

ði  1Þp ti ¼ T=2 1  cos ; m1

i ¼ 1; 2; . . . ; m

ð42Þ

where T is the time span. 4. DQ approximation of time derivatives For the DQ solution of Eq. (34), consider m sample time points in the time domain t1 = 0, t2, t3, . . ., tm = T, where T is the time domain of interest (i.e., time span). First note that

fyðtÞg ¼ ½ y1 ðtÞ y2 ðtÞ . . . y2n ðtÞ T ¼ ½ w1 ðtÞ w2 ðtÞ . . . wn ðtÞ u1 ðtÞ u2 ðtÞ . . . un ðtÞ T

ð43Þ

2501

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

From Eqs. (37) and (38), the DQ analogues of the first- and second-order time derivatives can be expressed as:

y_ k ðti Þ ¼

m X

€k ðt i Þ ¼ Aij yk ðt j Þ; y

j¼1

m X

Bij yk ðt j Þ;

k ¼ 1; 2; . . . ; 2n

ð44Þ

j¼1

Satisfying Eq. (34) at any sample time point t = ti, yields to the following:

9 8P m > > > > A y ðt Þ ij j 1 > > > > > > j¼1 > > > > > > > m > P > > > = < Aij y2 ðt j Þ >

9 9 8 8 y1 ðt1 Þ > y_ 1 ðt i Þ > > > > > > > > > > > > > = = < y2 ðt1 Þ > < y_ 2 ðt i Þ > j¼1 _ ¼ Ai1 . ¼ fyðt i Þg ¼ . > > > > > > .. . > > > > > > > > > > .. > > > > > > >. > ; ; > > : : . > > > > ðt Þ y y_ 2n ðt i Þ > > 1 2n m > >P > > > ; : Aij y2n ðtj Þ > j¼1

9 9 8 8 y1 ðt 2 Þ > y1 ðtm Þ > > > > > > > > > > > > > = = < y2 ðt 2 Þ > < y2 ðtm Þ > þ Ai2 . þ    þ Aim . > > > > .. .. > > > > > > > > > > > > ; ; : : y2n ðt 2 Þ y2n ðt m Þ

€ðt i Þg þ ½Cðti Þfyðt _ i Þg þ ½Kðti Þfyðt i Þg ¼ fFðti Þg; ½Mðt i Þfy i ¼ 1; 2; . . . ; m

ð45Þ

Using the quadrature rule, Eq. (44), the first-order derivative of the vector {y(t)} can be written as:

ð46Þ

-0.2

-0.2

0

0

n=5 nT=5, m=15

0.2

n=1 n=3 n=5 n=7

0.4

W cd / W cs

Wcd / W cs

0.2

v = 15.6 m/s

0.6

0.4

0.8

1

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.2

1

0

0.1

0.2

0.3

0.4

0

0 nT=1, m=18

0.7

0.8

0.9

1

n=1 n=2 n=4 n=6

0.4

v = 62.4 m/s

0.6 0.8

1.2 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.4

1

nT=6, m=7

0.8

1.2 0.2

v = 62.4 m/s

nT=4, m=9

0.6

1

0.1

nT=2, m=12

0.4

1

0

n=5 nT=1, m=18

0.2

Wcd / W cs

0.2

Wcd / W cs

0.6

-0.2

-0.2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t/T

t/T -0.2

-0.2

0

0 nT=1, m=12 n=1 n=2 n=3 n=5

0.4 0.6

1

0.6

1 1.2

1.4

1.4

1.6

1.6 0.1

0.2

0.3

0.4

0.5

t/T

0.6

0.7

0.8

0.9

1

v = 124.8 m/s

0.8

1.2

0

nT=1, m=12 nT=2, m=8 nT=3, m=7 nT=5, m=6

0.4

v = 124.8 m/s

0.8

1.8

n=5

0.2

Wcd / W cs

0.2

Wcd / W cs

0.5

t/T

t/T

1.4

v = 15.6 m/s

0.6

0.8

1.2

nT=1, m=51 nT=2, m=26 nT=4, m=16 nT=6, m=13

1.8

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

t/T

Fig. 2. Convergence of the mid-span displacements of a simply supported beam subjected to a moving force with respect to n, nT, and m.

1

2502

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

1.2

Percent error in W ( L / 2, 0.75 T )

Percent error in W ( L / 2, 0.5 T )

0 -0.2 -0.4 DQM

-0.6

n=1 n=3 n=5 n=7 n=9

-0.8 -1 -1.2 -1.4

1

24

25

26

27

28

29

30

31

0.4 0.2 0 -0.2 -0.4 15

32

Number of sample time points

Percent error in W ( L / 2, 0.75 T )

Percent error in W ( L / 2, 0.5 T )

0.2

18

19

20

21

22

23

24

0 -0.2 Wilson

-0.4

n=1 n=3 n=5 n=7 n=9

-0.6 -0.8 -1 -1.2 250

300

350

400

450

1.4 1.2 Wilson

1

n=1 n=3 n=5 n=7 n=9

0.8 0.6 0.4 0.2

0 100 120 140 160 180 200 220 240 260 280 300

500

Number of sample time points

0.2

Percent error in W ( L / 2, 0.75 T )

1.4

0 -0.2 -0.4 Newmark

-0.6

n=1 n=3 n=5 n=7 n=9

-0.8 -1 -1.2 -1.4 250

300

350

400

450

1.2 1

Newmark

0.8

n=1 n=3 n=5 n=7 n=9

0.6 0.4 0.2 0

-0.2 100 120 140 160 180 200 220 240 260 280 300

500

Number of sample time points

Number of sample time points Fig. 3. Convergence and error analysis of the DQ solution for mid-span displacement of a simply supported beam subjected to a moving force at the time level t = 0.5T for v = 15.6 m/s and comparison with the other two schemes (nT = 1, Dt = T/ (m  1)).

By similar calculations, it can be obtained that:

9 9 9 8 8 8 y1 ðt 1 Þ > y1 ðt 2 Þ > y1 ðt m Þ > > > > > > > > > > > > > > > > > > > = = = < y2 ðt 1 Þ > < y2 ðt 2 Þ > < y2 ðt m Þ > €ðt i Þg ¼ Bi1 . þ Bi2 . þ    þ Bim . fy > > > > .. > > .. .. > > > > > > > > > > > > > > > > > > ; ; ; : : : y2n ðt 1 Þ y2n ðt 2 Þ y2n ðtm Þ ð47Þ

Fig. 4. Convergence and error analysis of the DQ solution for mid-span displacement of a simply supported beam subjected to a moving force at the time level t = 0.75T for v = 50 m/s and comparison with the other two schemes (nT = 1, Dt = T/ (m  1)).

_ i Þg ¼ fyðt

m X

Aij fyðt j Þg;

€ðt i Þg ¼ fy

j¼1

m X

Bij fyðt j Þg

ð48Þ

j¼1

Substituting Eq. (48) into Eq. (45) yields to:

½Mðti Þ

m X

Bij fyðt j Þg þ ½Cðt i Þ

j¼1

¼ fFðt i Þg;

m X

Aij fyðt j Þg þ ½Kðti Þfyðt i Þg

j¼1

i ¼ 1; 2; . . . ; m

Eqs. (46) and (47) can be written in the following series forms as:

25

1.6

Number of sample time points

Percent error in W ( L / 2, 0.5 T )

17

1.8

0.4

-1.6 200

16

Number of sample time points

0.6

-1.4 200

n=1 n=3 n=5 n=7 n=9

0.6

-1.6 23

DQM

0.8

Eq. (49) can also be rewritten in the matrix form as:

ð49Þ

2503

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

2

e 11  ½K 6 e 21  6 ½K 6 6 . 6 . 4 .

e 12  ½K e 22  ½K .. .

... ... ...

e m1  ½ K e m2  . . . ½K

3 e 1m  8 fyðt 1 Þg 9 8 fFðt 1 Þg 9 ½K > > > > > > > > > > > 7> e 2m  7> = > = < fyðt 2 Þg > < fFðt 2 Þg > ½K 7 ¼ 7 . . .. 7> > > > . > > > .. > > > >. . 5> > > > ; > ; : : e mm  fyðt m Þg fFðt m Þg ½K

ð50Þ

5. Implementation of initial conditions and solution of resulting system of algebraic equations

e ij  are given by: where the 2n  2n sub-matrices ½ K

e ij  ¼ ½K



Bij ½Mðt i Þ þ Aij ½Cðt i Þ

i–j

ð51Þ

Bij ½Mðt i Þ þ Aij ½Cðt i Þ þ ½Kðt i Þ i ¼ j

Eq. (50) can also be written in the compact form as follows:

e f Y e g ¼ fe ½K Fg

Percent error in W ( L / 2, T )

Let the initial conditions of the problem are given by:

fyðt 1 Þg ¼ fyð0Þg ¼ fy0 g _ _ 1 Þg ¼ fyð0Þg ¼ fy_ 0 g fyðt

ð52Þ

ð54Þ

From Eq. (48), the derivative initial condition can be analogized

_ 1 Þg ¼ fy_ 0 g ¼ fyðt

2

m X

A1j fyðt j Þg

ð55Þ

j¼1

1.5

Substituting Eqs. (53) and (55) into Eq. (50), leads to the following system of algebraic equations:

DQM n=1 n=3 n=5 n=7 n=9

1 0.5

"

0 -0.5 -1

43

44

45

46

47

48

49

50

51

52

53

Number of sample time points 4

e ii  ½ K e id  ½K e di  ½ K e dd  ½K

#(

e ig fY e dg fY

)

( ¼

fe F ig f Fe d g

) ð56Þ

In this matrix equation, the subscripts i and d indicate the sample time points used for writing the quadrature analog of the initial conditions and the ordinary differential equations, respectively. By e i g; the matrix Eq. (56) is reduced eliminating the column vector f Y to the following system of algebraic equations:

b f Y e dg ¼ fb ½K Fg

ð57Þ

where

3.5

b  ¼ ½K e dd   ½ K e di ½ K e ii 1 ½ K e id  ½K 1 e di ½ K e ii  f e f Fb g ¼ f Fe d g  ½ K F ig

3 2.5 2

ð58Þ ð59Þ

Wilson

1.5

n=1 n=3 n=5 n=7 n=9

1 0.5 0

Table 1 Impact factors for a simply supported beam subjected to a moving concentrated force, obtained by the present mixed methodology and comparisons with those of conventional time step methods (n = 5, nT = 1). Note that the time steps for all schemes are the samea.

-0.5 -1 -1.5 -2

500 510 520 530 540 550 560 570 580 590 600

Number of sample time points 3.5 3

Percent error in W ( L / 2, T )

ð53Þ

as:

2.5

Percent error in W ( L / 2, T )

e  is the 2 nm  2 nm matrix and f Y e g and f e where ½ K F g are the 2 nm  1 vectors. After applying the initial conditions to Eq. (52), one can solve the resultant algebraic equations for unknown Ritz coefficients.

a

2.5 2 1.5 1

Tf/T

Velocity

Ref. [37]

Exact [51]

DQM

Wilson

Newmark

m

0.125 0.25 0.5 0.75 1 1.25 2

15.6 31.2 62.4 93.6 124.8 156 250

1.060 1.120 1.246 1.566 1.711 1.739 1.546

1.0604 1.1215 1.2585 1.5742 1.7057 1.7315 1.5462

1.0601 1.1213 1.2561 1.5734 1.7062 1.7118 1.5497

1.0281 1.1256 1.2847 1.5833 1.6570 1.6196 1.4827

1.0639 1.1262 1.2662 1.5635 1.6556 1.6013 1.4830

50 41 25 15 12 8 8

That is DtDQ = DtWilson = DtNewmark = T/(m  1), T = L/v.

Newmark n=1 n=3 n=5 n=7 n=9

0.5 0

Table 2 Impact factors for a simply supported beam subjected to a moving concentrated force obtained using the present mixed methodology and comparisons with those of other time step methods (n = 7, nT = 1, m = 100). Note that the time steps for all schemes are the same.a

-0.5 -1 -1.5 400 410 420 430 440 450 460 470 480 490 500

Number of sample time points Fig. 5. Convergence and error analysis of the DQ solution for mid-span displacement of a simply supported beam subjected to a moving force at the time level t = T for v = 125 m/s and comparison with the other two schemes (nT = 1, Dt = T/(m  1)).

a

Tf/T

Velocity

Exact [51]

DQM

Wilson

Newmark

0.125 0.25 0.5 0.75 1 1.25 2

15.6 31.2 62.4 93.6 124.8 156 250

1.0604 1.1215 1.2585 1.5742 1.7057 1.7315 1.5462

1.0603 1.1214 1.2587 1.5741 1.7057 1.7315 1.5463

1.0550 1.1237 1.2601 1.5749 1.7042 1.7336 1.5477

1.0615 1.1221 1.2578 1.5756 1.7034 1.7328 1.5470

That is DTDQ = DTWilson = DTNewmark = T/(m  1), T = L/v.

2504

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

Eq. (57) is a system of algebraic equations which can be solved using various direct or iterative methods. Since, the size of resulting algebraic equations is large (2 (nm  1)  2 (nm  1)), the direct methods (such as the Gaussian elimination method) are not efficient. Thus, the use of iterative methods (such as the Gauss–Seidel method) is recommended for the solution of the resultant algebraic equations.

of the beam are assumed to be clamped. Since the information regarding the forced vibration responses of a FG beam due to moving masses is not available, the results for homogeneous simply supported beams subjected to moving masses are first presented and compared with the existing results to confirm the reliability and accuracy of the present mixed method for the solution of moving mass problem. After that, the dynamic behavior of a FG beam due to a moving mass is investigated and the effects of some parameters such as velocity of the moving load, inertia of the moving load and ratio of material properties of the FG beam are examined. To solve Eq. (34) using the scheme described in Section 6, the time domain of interest is divided into nT equal DQ time element with m sample time points (per DQ time element). The total number of sample time points and the average time step can be obtained respectively as:

6. A step-by-step DQ in time When very long-term solutions are required or the size of resulting matrices in Eq. (57) is too large, it is more convenient to apply the DQM as a step-by-step time integration scheme to advance the solutions progressively over the time domain of interest. In this technique, the time domain of interest is first divided into several time elements. The DQM in then applied to each time element independently. Note that the solutions at the end of each DQ time element is used as the initial conditions for the next time element (for more details of this technique see Refs. [38–40]). This technique reduces considerably the computational cost, since DQM is applied to each time element independently and a smaller system of algebraic equations should be solved at each step. Also, note that each step in this technique consists of several sub-steps.

Mtot ¼ nT ðm  1Þ þ 1

ð60Þ

Dtav ¼ T=ðM tot  1Þ ¼ T=ðnT ðm  1ÞÞ

ð61Þ

where T is the time span. The numerical examples are also solved using classical finite difference methods such as the Newmark, and Wilson h methods. Parameters a = 0.25 and d = 0.5 are taken in the Newmark scheme and h = 1.4 in the Wilson h method [50]. 7.1. Vibration of a simply supported homogeneous beam due to a moving force

7. Numerical results and discussion A number of numerical examples are now presented to confirm the reliability, convergence, efficiency, and accuracy of the proposed mixed methodology. The first example is of moving force problem of a homogeneous simply supported beam. This problem has an analytical solution and the accuracy of the present mixed method is verified by comparing the calculated results with those obtained by analytical ones. In the second example, the dynamic responses of a simply supported FG beam subjected to a moving force are determined and compared with those of existing literature. The same problem is solved in the third example, but the end conditions

Consider a simply supported homogeneous beam with a span L of 10.16 cm, width b = 0.635 cm, thickness h = 0.635 cm, modulus of elasticity E = 2.068  1011 Pa, mass density q = 10686.9 kg/m3, moment of area I = 135.4832 mm4 and fundamental period Tf = 8.149  10–4 s, subjected to a moving force P = 4.45 N [36,37]. The present problem has an analytical solution [51] and the accuracy of the present mixed method is tested by comparing the calculated results with those of analytical solutions. The Rayleigh–Ritz method with n parameters and the DQM with nT number of time elements

Table 3 Maximum normalized deflections at the beam centre for different number of DQM time elements, nT, and sample time points (per DQM time elements) (n = 7). Method

nT

m

k = 0.2 (v = 222 m/s)

k = 0.5 (v = 198 m/s)

k=1 (v = 179 m/s)

k=2 (v = 164 m/s)

Pure alumina (v = 252 m/s)

Pure steel (v = 132 m/s)

Exp. lawa (v = 180 m/s)

Present

1 2 3 6 12 30 40 100

15 9 8 6 5 5 4 4

1.0333 1.0332 1.0332 1.0331 1.0340 1.0331 1.0333 1.0333 1.0344

1.142 1.1427 1.1427 1.1427 1.1436 1.1427 1.1429 1.1429 1.1444

1.2486 1.2483 1.2484 1.2483 1.2494 1.2484 1.2485 1.2486 1.2503

1.3360 1.3358 1.3358 1.3357 1.3368 1.3357 1.3360 1.3359 1.3376

0.9317 0.9316 0.9316 0.9315 0.9323 0.9315 0.9317 0.9317 0.9328

1.7302 1.7298 1.7299 1.7298 1.7313 1.7298 1.7301 1.7301 1.7324

1.2737 1.2736 1.2736 1.2734 1.2745 1.2735 1.2737 1.2737 1.2754

Ref. [22] a

Exponential law.

Table 4 Maximum normalized deflections at the beam centre obtained by various schemes (n = 7).

a b

Method

k = 0.2 (v = 222 m/s)

k = 0.5 (v = 198 m/s)

k=1 (v = 179 m/s)

k=2 (v = 164 m/s)

Pure alumina v = 252 m/s)

Pure steel (v = 132 m/s)

Exp. law (v = 180 m/s)

DQMa (nT = 1, m = 15) Wilsona (m = 15) Newmarka (m = 15) DQMb (nT = 1, m = 35) Wilsonb (m = 400) Newmarkb (m = 400) Ref. [22]

1.0333 1.0231 1.0216 1.0338 1.0339 1.0337 1.0344

1.1429 1.1315 1.1300 1.1434 1.1435 1.1434 1.1444

1.2486 1.2361 1.2345 1.2492 1.2493 1.2492 1.2503

1.3360 1.3229 1.3208 1.3365 1.3367 1.3365 1.3376

0.9317 0.9225 0.9212 0.9321 0.9322 0.9321 0.9328

1.7302 1.7130 1.7106 1.7310 1.7311 1.7309 1.7324

1.2737 1.2612 1.2593 1.2743 1.2744 1.2743 1.2754

DTDQ = DTWilson = DTNewmark = T/(m  1) = T/14 = 0.07143T, T = L/v. DTDQ = T/34 = 0.0294T, DTWilson = DTNewmark = T/399 = 0.0025T.

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

0.78 0.77

W ( L / 2, 0.5 T ) / W cs

0.76 0.75 v = 10 m/s, k=2

0.74

n=1 n=3 n=5 n=7 n=9

0.73 0.72 0.71 0.7 0.69 3

6

9

12

15

Number of sample time points 0.605 0.6

W ( L / 2, 0.5 T ) / W cs

0.595 0.59 0.585 v = 1 m/s, k=0.2

0.58

n=1 n=3 n=5 n=7 n=9

0.575 0.57 0.565 0.56 0.555 3

6

9

12

15

Number of sample time points 0.76

W ( L / 2, 0.5 T ) / W cs

0.74 v = 100 m/s, Exp. law

0.72

n=1 n=3 n=5 n=7 n=9

0.7 0.68 0.66 0.64 10

11

12

13

14

15

16

17

18

19

7.2. Vibration of a simply supported FG beam due to a moving force

Fig. 6. Convergence of the DQ solution for mid-span displacement of a clamped– clamped beam subjected to a moving force at the time level t = 0.5T for different values of v, k, and n (nT = 1, Dt = T/(m  1)).

with m sample time points (per each DQ time element) are used in the mixed method for the solution of the present problem. The Rayleigh–Ritz trial functions are assumed to be the normal modes (mode shapes) of the simply supported beam, i.e.,

jp x L

for different values of n, nT, and m. The use of Lagrange interpolations in the time direction enables us to find a continuous representation of the approximate solutions, as it can be seen from the results shown in Fig. 2. The converging trend of solutions with increasing n, nT, and m is obvious in Fig. 2. It can be seen that the use of only 3–5 Ritz test functions are sufficient to achieve accurate solutions. Besides, by increasing the number of time elements (nT), a smaller number of sample time points per DQ time element are required to obtain converged solutions. Since the employed DQM is a Lagrange-polynomial-based DQM, thus, it is the shape of dynamic responses which determines the convenient values of nT and m for insuring the stability and accuracy of solutions. For example, the smoother shape of dynamic responses; then the smaller values of nT and m. To demonstrate the rate of convergence of the present mixed methodology, the percent error in mid-span displacements (defined as (wNumeric  wExact)/wExact  100) is calculated for different values of n, m and v. In Figs. 3–5, the percent error in solutions is plotted against the number of sample time points for different values of n and v at the time levels t = 0.5T, t = 0.75T, and t = T, respectively (where T is the traveling time of the force moving from the left end to the right end of the beam). The DQ solutions are obtained using a single time element (i.e., nT = 1). The solutions of the Newmark and Wilson methods are also shown for comparison. It can be observed that the DQ method produces much better accuracy than the Newmark and Wilson methods using a considerably much smaller number of sample time points and therefore requires relatively little computational effort. It can also be observed that all the time integration schemes produce accurate solutions by using n = 3 only. Also, note that the solutions of all schemes with n = 1 are not acceptable for accuracy. The results for the impact factor, defined as the ratio of maximum dynamic deflection to maximum static deflection at the beam centre, are tabulated in Tables 1 and 2. The results are given for two situations, i.e., when a rather large time step sizes are employed (see Table 1) and when a rather small time step sizes are used (see Table 2). The results of the present mixed method are obtained using a single DQM time element. The results of existing literature and those of single-step methods are also cited for comparison. It can be seen that for both situations, the DQM gives better accuracy than the Newmark and Wilson h methods. This is actually due to the fact that the DQM is a high-order method, while the finite difference schemes such as the Wilson and Newmark methods are second-order methods only. It can also be observed that smaller number of sampling time points (i.e., larger time step sizes) are allowed to be used when using the DQ time integration scheme.

20

Number of sample time points

/j ðxÞ ¼ sin

2505

ð62Þ

A convergence study is first performed to show the rapid convergence and accuracy of the mixed scheme. The dynamic responses at the beam centre, wcd, are evaluated for different values of v (velocity of the moving point load) and normalized by the static deflection wcs = PL3/48EI of the beam under a point load P at midspan. Fig. 2 demonstrates the convergence behavior of solutions

Consider a simply supported FG beam with a span L of 20 m, width b = 0.4 m, thickness h = 0.9 m, subjected to a moving point load P = 100 kN [22]. Assume that the FG beam is composed of steel and alumina (Al2O3) and its properties varies through the thickness of the beam according to power-law and exponential law (see Eqs. (1) and (2)). Thus, the bottom surface of the beam is pure steel, while the top surface of the beam is pure alumina. The material properties of steel and alumina are given by:

Alumina : Steel :

Et ¼ 390 GPa; qt ¼ 3960 kg=m3

Eb ¼ 210 GPa; qb ¼ 7800 kg=m3

For the present test problem, the Rayleigh–Ritz trial functions are assumed to be:

/j ðxÞ ¼ sin

jp jp x; wj ðxÞ ¼ cos x L L

ð63Þ

2506

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

The dynamic responses at the beam centre are evaluated for various travel speeds (v) and material distributions and normalized by the static deflection wcs = PL3/48EbI of the fully steel beam under a point load p at mid-span. Table 3 gives the maximum normalized deflections at the beam centre for different values of v (velocity of moving load) and k (power-law exponent). The results of Ref. [22] are also included for comparison purposes. Excellent agreements are achieved between the results of the present analysis and those of Ref. [22]. Note that the results of Ref. [22] were obtained using 14 terms of polynomial functions series and 500 time steps. From the results shown in Table 3, it can also be observed that the solutions exhibit uniform convergence without instability for an increase in the number of DQM time elements (nT) and number of sampling time points (m). The resultant dynamic equations of the present problem are also solved using the Newmark and Wilson h methods. In Table 4, the results of DQ time integration scheme are compared with those of single time step methods. It can be seen that the DQM gives much better accuracy than the Newmark and Wilson methods using a much smaller number of sample time points (i.e., using a much larger time step sizes). Note that the DQ results with only m = 35 are comparable in accuracy to those of single-step methods with m = 400. 7.3. Vibration of a clamped–clamped FG beam due to a moving force Consider an FG beam with clamped end conditions subjected to a moving force. In this case, the Rayleigh–Ritz trial functions are assumed to be in the following forms [48]:

/j ðxÞ ¼ ðx=LÞiþ1  2ðx=LÞiþ2 þ ðx=LÞiþ3 ;

wj ðxÞ ¼ sin

jp x L

ð64Þ

It is not difficult to verify that the above trial functions satisfy the homogeneous form of boundary conditions of the clamped– clamped beam. The material properties of the beam are assumed to be the same as those employed in Section 7.2. Similar to Section 7.2, the beam centre deflections are evaluated for various values of v and k and normalized by the static deflection wcs = PL3/192EbI of the fully steel beam under a point load p at mid-span. Fig. 6 illustrates the results at the time level t = 0.5T for different values of n and m. A single DQM time element is used in the time domain to obtain the solutions. The converging trend of solutions with increasing n and m is obvious in Fig. 6. By comparing the results with those shown in Figs. 3–5, it is concluded that the rate of convergence is slower with clamped–clamped boundary conditions. In other words, the number of Ritz terms used for clamped–clamped boundary conditions should be greater then the one used for simply supported boundary conditions. It may

be seen from Fig. 6 that the results of good accuracy for clamped–clamped beam can be obtained if n P 5. The results of the impact factor are given in Table 5. From the results shown in Table 5, it can be observed that the DQ solutions with only 17 sample time points are comparable in accuracy to those of other time step methods (i.e., Newmark and Wilson methods) with 200 sample time points. These results confirm the highefficiency and accuracy of the DQM time integration scheme again. 7.4. Vibration of a simply supported homogeneous beam due to a moving mass Since there is no information regarding the forced vibration of FG beams subjected to moving masses, this numerical example is presented to validate the proposed formulation for solving the moving mass problem by comparing the numerical results with those of existing literature. Consider a simply supported homogeneous beam subjected to a moving mass. The material properties of the beam are assumed to be the same as those employed in Section 7.1. The rate of convergence of the impact factor with different number of time elements, traveling speeds and point mass values is used to demonstrate the effectiveness of the proposed mixed methodology. In solving the present problem, it is found that the DQ time integration scheme may encounter some convergence problem and numerical instability if m (number of sample time points per DQM time element) be rather large. This undesirable feature has also been reported by some researchers in solving static problems (see Ref. [26] and references therein). In solving the second-order initial-value problems using DQM, Fung [34] also addressed this difficulty. He showed that the DQ time integration scheme is not unconditionally stable in general and to achieve unconditionally stable algorithms using the DQM, the number of sample time points should be rather small. Thus, the numerical results presented in this section are obtained using small m-values. Fig. 7 shows the convergence of the mid-span deflections of a simply supported beam subjected to a moving mass for three different point mass values. The mass values are normalized with respect to the beam mass and presented in terms of a mass ratio Mratio = Mmoving/Mbeam. From the results shown in Fig. 7, it can be observed that as the mass ratio increases a larger number of DQ time elements are required for convergence and to achieve accurate solutions. It can also be seen that the Rayleigh–Ritz method with only 2–5 terms produces satisfactory solutions. Table 6 presents the normalized deflections at the beam centre at the time level t = 0.5T for different number of Ritz terms and point mass values. The symbols Tf and T in Table 6, respectively, denote the fundamental period of the beam and the traveling time of the moving mass. An excellent convergence trend of solutions can

Table 5 Maximum normalized deflections at the beam centre for different number of DQM time elements, nT, and sample time points (per DQM time elements (n = 7). Method

nT

m

k = 0.2 (v = 222 m/s)

k = 0.5 (v = 198 m/s)

k=1 (v = 179 m/s)

k=2 (v = 164 m/s)

Pure alumina (v = 252 m/s)

Pure steel (v = 132 m/s)

Exp. law (v = 180 m/s)

DQMa

1 2 3 4 5 35 – – – –

17 13 9 8 7 4 17 17 200 200

0.8283 0.8281 0.8290 0.8287 0.8284 0.8281 0.8333 0.8224 0.8285 0.8282

0.9168 0.9166 0.9176 0.9172 0.9169 0.9166 0.9220 0.9100 0.9170 0.9166

1.0022 1.0018 1.0030 1.0026 1.0022 1.0019 1.0074 0.9945 1.0024 1.0020

1.0708 1.0706 1.0717 1.0713 1.0709 1.0705 1.0773 1.0632 1.0711 1.0706

0.7468 0.7467 0.7475 0.7472 0.7469 0.7466 0.7513 0.7415 0.7470 0.7467

1.3886 1.3880 1.3897 1.3890 1.3885 1.3881 1.3959 1.3780 1.3888 1.3882

1.0202 1.0202 1.0212 1.0208 1.0205 1.0200 1.0269 1.0133 1.0206 1.0202

Wilsona Newmarka Wilsonb Newmarkb a b

DTDQ = DTWilson = DTNewmark = T/(m  1) = T/16 = 0.0625T, T = L/v. DTWilson = DTNewmark = T/199 = 0.0050T.

2507

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

be observed. The results show that the normalized deflections converge to their final values with n = 7 (up to three decimal points) regardless of point mass values. Note that the rate of convergence is not affected by the point mass values. Table 7 presents the impact factor at the beam centre for three different normalized time periods Tf/T, in relation to seven different DQM time discretizations. The results of Ref. [52] are also included for comparison. An excellent agreement is observed. It can also be seen that as the mass ratio increases, the impact factor convergence requires more DQM time elements. The results shown in Tables 6 and 7 are also suggested that the rate of convergence is more sensitive

to the number of DQM time elements than the number of Ritz terms. 7.5. Vibration of a simply supported FG beam due to a moving mass In this section, the dynamic behavior of a simply supported FG beam subjected to a moving mass is analyzed. In the analysis performed, the effects of various parameters affecting the dynamic behavior of the FG beam such as velocity and inertia of the moving load and material properties of the beam are investigated. For convenience of comparison, the dimensions and material properties of

-0.2

-0.2

0

0 nT = 40, m = 4 n=1 n=3 n=5 n=7

W cd / W cs

0.4

Mratio = 0.5

0.8

0.6

1

1.2

1.2

1.4

1.4

1.6 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.6 0

1

0.1

0.2

0.3

0.4

t/T

0.6

0.7

0.8

0.9

1

0.6

0.7

0.8

0.9

1

0.7

0.8

-1 nT = 40, m=4

-0.5

n = 7, m=4

-0.5

n=1 n=2 n=3 n=5

Mratio = 1

nT = 15 nT = 25 nT = 35 nT = 50

0

W cd / W cs

0

W cd / W cs

0.5

t/T

-1

0.5

Mratio = 1

0.5

1

1

1.5

1.5

2

2 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

t/T

0.5

t/T

-2

-2.5 -2

-1.5

n=7, m=4

nT =65, m=4

-1.5

n=1 n=2 n=3 n=5

-0.5

nT = 25 nT = 45 nT = 65 nT = 85

-1

Mratio = 2

W cd / W cs

-1

W cd / W cs

Mratio = 0.5

0.8

1

0.1

nT = 10 nT = 20 nT = 30 nT = 40

0.4

0.6

0

n = 7, m = 4

0.2

W cd / W cs

0.2

0 0.5

-0.5

Mratio = 2

0 0.5 1

1 1.5 1.5

2

2 0

0.1

0.2

0.3

0.4

0.5

t/T

0.6

0.7

0.8

0.9

1

2.5 0

0.1

0.2

0.3

0.4

0.5

0.6

0.9

1

t/T

Fig. 7. Convergence of the mid-span displacements of a simply supported beam subjected to a moving mass force with respect to n and nT. (v = 62.4 m/s).

2508

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

Table 6 Maximum normalized deflections at the beam centre at the time level t = 0.5T for different number of Ritz terms and point mass values (nT = 50, m = 4). Mratio

Tf /T

n=3

n=5

n=7

n=9

n = 11

n = 13

n = 15

n = 17

0.5

0.5 1.0 1.5

1.4143 1.0120 0.5680

1.4160 1.0143 0.5717

1.4166 1.0146 0.5722

1.4168 1.0148 0.5725

1.4168 1.0149 0.5726

1.4169 1.0150 0.5727

1.4169 1.0150 0.5726

1.4169 1.0150 0.5726

1.5

0.5 1.0 1.5

1.4097 0.6919 0.3624

1.4108 0.6966 0.3676

1.4108 0.6975 0.3688

1.4109 0.6978 0.3692

1.4109 0.6979 0.3694

1.4109 0.6979 0.3694

1.4109 0.6979 0.3694

1.4109 0.6979 0.3694

3.0

0.5 1.0 1.5

1.1862 0.4729 0.2359

1.1893 0.4780 0.2401

1.1896 0.4793 0.2413

1.1896 0.4797 0.2417

1.1896 0.4799 0.2419

1.1896 0.4800 0.2420

1.1895 0.4800 0.2420

1.1895 0.4800 0.2420

Table 7 Impact factors at the beam centre for different number of DQM time elements, nT, and point mass values (n = 7, m = 4). Mratio

Tf/T

nT = 25

nT = 50

nT = 75

nT = 100

nT = 125

nT = 150

nT = 200

Ref. [52]

0.5

0.5 1.0 1.5

1.4084 2.1141 2.2957

1.4157 2.0258 2.2589

1.4171 2.0373 2.2657

1.4176 2.0455 2.2667

1.4180 2.0482 2.2677

1.4183 2.0489 2.2682

1.4186 2.0492 2.2689

1.418 2.048 2.271

1.0

0.5 1.0 1.5

1.5639 2.6647 2.6653

1.5684 2.4794 2.7441

1.5706 2.4612 2.6729

1.5716 2.4567 2.7046

1.5723 2.4558 2.7023

1.5727 2.4554 2.7007

1.5730 2.4554 2.6990

1.571 2.453 2.690

2.0

0.5 1.0 1.5

1.9659 3.3652 2.7541

1.9041 3.4063 3.4427

1.8845 3.9793 3.3687

1.8902 3.4028 3.7917

1.8795 3.4072 3.6099

1.8752 3.4057 3.6559

1.8709 3.4095 3.7097

1.857 3.404 3.709

Table 8 Maximum normalized deflections at the beam centre by moving force model.

v (m/s)

Pure alumina

k = 0.2

k = 0.5

k=1

k=2

k=5

k = 10

k = 20

Pure steel

Exp. law

20 40 60 80 100 125 150 175 200 225 250 275 300

0.5656 0.5884 0.6308 0.5824 0.6866 0.7852 0.8540 0.8980 0.9209 0.9304 0.9321 0.9262 0.9125

0.6252 0.6476 0.6872 0.7119 0.8260 0.9259 0.9890 1.0206 1.0325 1.0330 1.0227 1.0035 0.9760

0.7048 0.7511 0.7175 0.8538 0.9732 1.0715 1.1223 1.1409 1.1428 1.1293 1.1042 1.0679 1.0331

0.7706 0.8411 0.8155 0.9941 1.1155 1.2058 1.2417 1.2495 1.2370 1.2077 1.1641 1.1222 1.0781

0.8167 0.9043 0.9290 1.1152 1.2354 1.3123 1.3356 1.3313 1.3022 1.2555 1.2061 1.1550 1.0971

0.8903 0.9523 1.0433 1.2336 1.3491 1.4104 1.4219 1.4014 1.3561 1.2975 1.2397 1.1751 1.1014

0.9466 0.9893 1.1332 1.3275 1.4395 1.4912 1.4949 1.4616 1.4038 1.3413 1.2754 1.1974 1.1177

1.0009 1.0252 1.2218 1.4208 1.5291 1.5736 1.5692 1.5241 1.4572 1.3888 1.3132 1.2247 1.1410

1.0993 1.0862 1.3911 1.5983 1.6989 1.7316 1.7096 1.6454 1.5640 1.4810 1.3825 1.2816 1.1906

0.7868 0.8568 0.8266 1.0090 1.1337 1.2274 1.2656 1.2747 1.2633 1.2345 1.1911 1.1486 1.1042

Table 9 Maximum normalized deflections at the beam centre by moving mass model.

v (m/s)

Pure alumina

k = 0.2

k = 0.5

k=1

k=2

k=5

k = 10

k = 20

Pure steel

Exp. law

20 40 60 80 100 125 150 175 200 225 250 275 300 325 350

0.5547 0.5718 0.5965 0.6269 0.7439 0.8525 0.9275 0.9848 1.0331 1.0738 1.0913 1.1077 1.1223 1.1182 1.0971

0.6305 0.6680 0.6134 0.7611 0.8863 0.9942 1.0684 1.1257 1.1651 1.1844 1.2019 1.2119 1.1946 1.1638 1.1229

0.6963 0.7579 0.7267 0.9070 1.0368 1.1427 1.2169 1.2648 1.2862 1.3040 1.3099 1.2847 1.2414 1.1894 1.1355

0.7568 0.8292 0.8570 1.0510 1.1798 1.2856 1.3524 1.3825 1.3999 1.4066 1.3776 1.3252 1.2635 1.2016 1.1750

0.8305 0.8723 0.9732 1.1740 1.3003 1.4020 1.4530 1.4736 1.4853 1.4618 1.4060 1.3366 1.2616 1.2277 1.1914

0.8937 0.8987 1.0901 1.2927 1.4173 1.5098 1.5435 1.5600 1.5492 1.4967 1.4208 1.3372 1.2858 1.2434 1.1968

0.9419 0.9194 1.1829 1.3882 1.5133 1.5983 1.6261 1.6369 1.6074 1.5360 1.4493 1.3642 1.3193 1.2710 1.2139

0.9880 0.9388 1.2749 1.4840 1.6100 1.6877 1.7131 1.7161 1.6683 1.5822 1.4846 1.4116 1.3605 1.3038 1.2390

1.0689 1.0799 1.4513 1.6681 1.7951 1.8597 1.8812 1.8620 1.7828 1.6751 1.5636 1.5057 1.4425 1.3682 1.2957

0.7701 0.8459 0.8702 1.0688 1.2015 1.3107 1.3815 1.4149 1.4334 1.4434 1.4161 1.3641 1.3015 1.2363 1.2112

the FG beam studied in this section are chosen to be identical to those of Ref. [22] (see Section 7.2). Note that the weight and mass

of the moving load are, respectively, Mg = 100 kN and M = 10.1937  103 kg.

2509

0.9

1

0.8

0.9

0.7

0.8 0.7

0.6

w (vt, t) / w cs

w (vt, t) / w cs

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

0.5 0.4 v = 20 m/s

0.3

Moving force (k = 0.5) Moving mass (k = 0.5) Moving force (k = 5) Moving mass (k = 5)

0.2

0.6 0.5 0.4

v = 35 m/s Moving force (k = 0.5) Moving mass (k = 0.5) Moving force (k = 5) Moving mass (k = 5)

0.3 0.2

0.1

0.1

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t/T

1.2

1.4

1

1.2 1

0.8

w (vt, t) / w cs

w (vt, t) / w cs

t/T

0.6 0.4

v = 60 m/s

0.6 v = 80 m/s

0.4

Moving force (k = 0.5) Moving mass (k = 0.5) Moving force (k = 5) Moving mass (k = 5)

0.2

0.8

Moving force (k = 0.5) Moving mass (k = 0.5) Moving force (k = 5) Moving mass (k = 5)

0.2

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

t/T

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t/T

Fig. 8. Normalized dynamic responses of the beam under moving load for different values of

To investigate the effect of inertia of the moving load, the maximum normalized deflections at the beam centre are evaluated for different values of v (velocity of the moving load) and k (power-law exponent) for both the moving force and the moving mass models. The results are given in Tables 8 and 9. In the results shown in Tables 8 and 9, the maximum magnitudes of maximum normalized deflections of the FG beam centre are bolded. Note that when k approaches to zero, the material properties of the FG beam approaches to those of pure alumina and when k approaches to infinity, the FG beam material properties approaches to those of pure steel. Therefore, it is expected that the deflections of the FG beam approach to those of pure alumina and steel beam as k approaches to zero and infinity, respectively. This phenomenon can be clearly seen from the results shown in Tables 8 and 9. On the other hand, the deflections of the FG beam can be controlled by choosing the proper values of k. By comparing the results of Table 8 (moving force model) and those of Table 9 (moving mass model), one can easily see the effects of inertia of moving load. As it can be seen, the moving force model may predict larger or smaller deflections than the moving mass model when the velocity of the moving load is low (i.e., when v < 60 m/s). This behavior is independent of the material properties of the FG beam as seen in Tables 8 and 9. Nevertheless, the moving mass model predict larger deflections than the moving force model at high speed (i.e., when v P 60 m/ s). In other words, while the effect of inertia of the moving load can be neglected at low speeds it cannot be neglected at high speeds regardless of the material properties of the FG beam. Another interesting and important results of the present study is that the critical velocity (the velocity at which the maximum magnitude of maximum deflection of the beam to be occurred) of the FG beam can also be controlled by choosing the suitable values

0.2

v (velocities of moving load) and k (power-law exponent).

of k. From the results shown in Tables 8 and 9, it can also be seen that the moving force model can predict lower critical speeds than the moving mass model for maximum normalized deflections. For example when k = 0.2, the critical velocity predicted by the moving force model is approximately vcr = 225 m/s which is much smaller than the corresponding one by the moving mass model (vcr = 275 m/s). Fig. 8 illustrates the dynamic responses of the FG beam under the moving load for different moving load speeds and material properties of the beam. It is seen that as the velocity of moving load increases, the shape of dynamic responses becomes smoother. It is also seen that the shape of dynamic responses obtained from the moving force model differs from that of the moving mass model. From the results shown in Fig. 8, one can also clearly see the effect of inertia of the moving load, i.e., as the velocity of moving load increases, the moving mass model predict larger maximum displacements than the moving force model. This effect can be neglected when velocity of moving load is low, as seen in Fig. 8. As pointed out earlier, this effect becomes more crucial as the velocity of the moving load increases. The axial displacements of the beam centre are also calculated for different values of v and k and normalized by the static displacement ucs = PL/AEb (where Eb is the Young modulus at the bottom surface of the beam). Tables 10 and 11 give the results for the moving force and moving mass models, respectively. It can be observed that the larger the velocities of the moving load, the larger the axial displacements of the FG beam. Besides, as k increases the axial displacements of the FG beam vary from the displacements of the pure alumina beam to those of the pure steel beam as expected. From the results shown in Tables 10 and 11, one can also see the effects of inertia of the moving load on the axial displacements of the FG

2510

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

Table 10 Maximum normalized axial displacements at the beam centre by moving force model.

v (m/s)

Pure alumina

k = 0.2

k = 0.5

k=1

k=2

k=5

k = 10

k = 20

Pure steel

Exp. law

15 30 75 150 250 300 350

0.5865 0.5975 0.6339 0.6887 0.8505 0.9382 1.0069

0.6261 0.6418 0.6899 0.7870 0.9697 1.0636 1.1471

0.6706 0.6886 0.7224 0.8137 1.1007 1.1760 1.5260

0.7131 0.7325 0.8018 0.8899 1.2156 1.4753 1.8455

0.7555 0.7806 0.8543 1.0022 1.3184 1.7903 2.1170

0.8240 0.8531 0.9531 1.1772 1.6281 2.1416 2.4285

0.8984 0.9214 1.0446 1.3013 1.9307 2.4364 2.6874

0.9688 1.0062 1.1347 1.4355 2.2189 2.7225 2.9299

1.1086 1.1545 1.2889 1.6981 2.7430 3.2326 3.3407

0.7252 0.7465 0.8126 0.8994 1.2337 1.4766 1.8603

Table 11 Maximum normalized axial displacements at the beam centre by moving mass model.

v (m/s)

Pure alumina

k = 0.2

k = 0.5

k=1

k=2

k=5

k = 10

k = 20

Pure metal

Exp. law

15 30 75 150 250 300 350

0.5889 0.6139 0.6793 1.0194 1.5966 1.7310 1.6570

0.6249 0.6521 0.7168 1.0987 1.7001 1.6027 1.6448

0.6784 0.6816 0.7235 1.0127 1.5842 1.5574 1.8288

0.7 0.7393 0.7610 1.3250 1.5686 1.6553 2.3132

0.7551 0.8020 0.8573 1.5747 1.4395 2.0543 2.8765

0.8187 0.8509 0.9579 1.8157 1.7767 2.5091 3.4269

0.9040 0.9333 1.0072 2.0055 2.1051 2.9021 4.0131

0.9825 1.0176 1.1019 2.1646 2.4487 3.3311 4.6729

1.1162 1.1681 1.4324 2.4087 3.0977 4.1892 5.7763

0.7311 0.7528 0.7713 1.3395 1.6359 1.6606 2.3860

beam. It can be seen that the inertia of the moving load has considerable influences on the longitudinal vibrations of the FG beam, especially when the velocity of the moving load is high. From the results of Tables 8–11, it can be observed that the influence of inertia of the moving load on the axial displacements of the FG beam is more crucial than that on the transverse displacements of the FG beam. Therefore, this effect should be considered in the analysis and design of FG beams subjected to moving loads.

in the design and analysis of FG structures subjected to moving dynamic loads. Another interesting result of the present study is that the critical velocity of the beam/FG beam predicted by moving force model is smaller than the one predicted by moving mass model. Besides, the critical velocity of the FG beam can be controlled by choosing the convenient values of k (power-law exponent). Appendix A

8. Conclusion In this paper, a mixed method is presented to analyze the dynamic behavior of FG beams subjected to moving loads. The efficiency, accuracy, and convergence of the proposed mixed method are investigated and analyzed in detail. The mixed method enjoys the simplicity of the Rayleigh–Ritz method and high-efficiency and accuracy of the DQM. It is shown that the proposed mixed methodology can predict the dynamic behavior of FG beams accurately. Numerical comparisons between the DQ time integration method and single time step method such as the Newmark and Wilson methods reveal that the DQM is a promising and effective tool for handling the linear systems of ordinary differential equations. Besides, the DQM can produce better accuracy than the time step methods using much larger time step sizes. The effects of velocity and inertia of the moving load, and material properties of the FG beam on the dynamics of the sys-

The axial force N, bending moment M, and transverse shear force Q are related to the normal strain e0 = u0,x and the flexural curvature kx = w0,xx by:



N M



 ¼

Axx

Bxx

Bxx

Dxx



e0 kx

 ;

Q ¼ M ;x ¼ Bxx e0;x  Dxx kx;x

ð65Þ

where Axx, Bxx and Dxx are given in Eq. (10) and (),x denotes the partial derivative with respect to x. For simply supported (SS) and clamped–clamped (CC) beams, the boundary conditions are respectively as:

N ¼ w0 ¼ M ¼ 0 at x ¼ 0; L u0 ¼ w0 ¼ w0;x ¼ 0 at x ¼ 0; L

ð66Þ ð67Þ

Each of the Ritz trial functions (/j(x) and wj(x) in Eqs. (14) and (15)) should satisfy at least geometric boundary conditions, i.e.,

/j ð0Þ ¼ /j ðLÞ ¼ 0 and wj;x ð0Þ ¼ wj;x ðLÞ ¼ 0 for SS beams

ð68Þ

/j ð0Þ ¼ /j;x ð0Þ ¼ /j ðLÞ ¼ /j;x ðLÞ ¼ 0 and wj ð0Þ ¼ wj ðLÞ ¼ 0 for CC beams

ð69Þ

tem are also investigated. It is found that the above-mentioned parameters have considerable effects on the dynamic behavior of the system. It is revealed that the effect of inertia of the moving load can be ignored at low speeds. Nevertheless, this effect can become more critical at high speeds. Moreover, the effect of inertia of the moving load on the axial vibrations of the beam is more crucial than the corresponding one on the transverse vibrations of the beam. Therefore, this effect should be considered

When the boundary conditions of the beam are simple (i.e., for SS beams), the normal modes of the beam can be easily chosen as the Ritz trial function (see example 1, Section 7.1). Otherwise (for CC beams, in which the natural frequencies are obtained by solving a nonlinear equation), the polynomial functions which satisfy the geometric boundary conditions of the problem can be chosen as the Ritz trial function (see example 3, Section 7.3).

S.M.R. Khalili et al. / Composite Structures 92 (2010) 2497–2511

References [1] Suresh S, Mortensen A. Fundamentals of functionally graded materials. London: IOM Communications Ltd.; 1998. [2] Suresh S, Mortensen A. Modeling and design of multi-layered and graded materials. Prog Mater Sci 1997;42:243–51. [3] Chakraborty A, Gopalakrishnan S, Reddy JN. A new beam finite element for the analysis of functionally graded materials. Int J Mech Sci 2003;45:519–39. [4] Sankar BV. An elasticity solution for functionally graded beams. Compos Sci Technol 2001;61:689–96. [5] Sankar BV, Taeng JT. Thermal stresses in functionally graded beams. AIAA J 2002;40:1228–32. [6] Venkataraman S, Sankar BV. Elasticity solution for stresses in a sandwich beam with functionally graded core. AIAA J 2003;41:2501–5. [7] Zhu H, Sankar BV. A combined Fourier series – Galerkin method for the analysis of functionally graded beams. ASME J Appl Mech 2004;71:421–3. [8] Shi ZF, Chen Y. Functionally graded piezoelectric cantilever beam under load. Arch Appl Mech 2004;74:237–47. [9] Nirmala K, Upadhyay PC, Prucz J, Loyns D. Thermoelastic stresses in composite beams with functionally grade layer. J Reinf Plast Compos 2005;24:1965–77. [10] Ching HK, Yen SC. Meslless local Petrov–Galerkin analysis for 2D functionally graded elastic solids under mechanical and thermal loads. Compos Part B – Eng 2005;36:223–40. [11] Ching HK, Yen SC. Transient thermoelastic deformations of 2D functionally graded beams under nonuniformly connective heat supply. Compos Struct 2006;73:381–93. [12] Lu CF, Chen WQ. Free vibration of orthotropic functionally graded beams with various end conditions. Struct Eng Mech 2005;13:1430–7. [13] Wu L, Wang QS, Elishakoff I. Semi-inverse method for axially functionally graded beams with an anti-symmetric vibration mode. J Sound Vib 2005;284:1190–202. [14] Aydogdu M, Taskin V. Free vibration analysis of functionally graded beams with simply supported edges. Mater Des 2007;28:1651–6. [15] Yang J, Chen Y. Free vibration and buckling analysis of functionally graded beams with edge cracks. Compos Struct 2008;83:48–60. [16] Kapuria S, Bhattacharyya M, Kumar AN. Bending and free vibration response of layered functionally graded beams: a theoretical model and its experimental validation. Compos Struct 2008;82:390–402. [17] Zhong Z, Yu T. Analytical solution of a cantilever functionally graded beam. Compos Sci Technol 2007;67:481–8. [18] Ying J, Lü CF, Chen WQ. Two-dimensional elasticity solutions for functionally graded beams resting on elastic foundations. Compos Struct 2008;84:209–19. [19] Li XF. A unified approach for analyzing static and dynamic behaviors of functionally graded Timoshenko and Euler–Bernoulli beams. J Sound Vib 2008;318:1210–29. [20] Rahmani O, Khalili SMR, Malekzadeh K, Hadavinia H. Free vibration analysis of sandwich structures with a flexible functionally graded syntactic core. Compos Struct 2009. doi:10.1016/j.compstruct.2009.05.007. [21] Yang J, Chen Y, Xiang Y, Jia XL. Free and forced vibration of cracked inhomogeneous beams under an axial force and a moving load. J Sound Vib 2008;312:166–81. [22] S ß imsßek M, Kocatürk T. Free and forced vibration of a functionally graded beam subjected to a concentrated moving harmonic load. Compos Struct 2009;90:465–73. [23] Bellman RE, Casti J. Differential quadrature and long term integrations. J Math Anal Appl 1971;34:235–8. [24] Bellman RE, Kashef BG, Casti J. Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations. Comput Phys 1972;10:40–52. [25] Shu C. Differential quadrature and its application in engineering. New York: Springer Publication; 2000. [26] Bert CW, Malik M. Differential quadrature method in computational mechanics: a review. ASME J Appl Mech Rev 1996;49:1–28. [27] Wu TY, Liu GR. The generalized differential quadrature rule for initial-value differential equations. J Sound Vib 2000;233(2):195–213.

2511

[28] Wu TY, Liu GR. Numerical solution for differential equations of duffing-type non-linearity using the generalized differential quadrature. J Sound Vib 2000;237:805–17. [29] Shu C, Yao KS. Block-marching in time with DQ discretization: an efficient method for time-dependent problems. Comput Methods Appl Mech Eng 2002;191:4587. 97. [30] Tanaka M, Chen W. Coupling dual reciprocity BEM and differential quadrature method for time-dependent diffusion problems. Appl Math Model 2001;25:257–68. [31] Tanaka M, Chen W. Dual reciprocity BEM applied to transient elastodynamic problems with differential quadrature method in time. Comput Methods Appl Mech Eng 2001;190:2331–47. [32] Fung TC. Solving initial value problems by differential quadrature method – part 1: first-order equations. Int J Numer Methods Eng 2001;50:1411–27. [33] Fung TC. Solving initial value problems by differential quadrature method – part 2: second- and higher-order equations. Int J Numer Methods Eng 2001;50:1429–54. [34] Fung TC. Stability and accuracy of differential quadrature method in solving dynamic problems. Comput Methods Appl Mech Eng 2002;191:1311–31. [35] Zhong H, Lan M. Solution of nonlinear initial-value problems by the splinebased differential quadrature method. J Sound Vib 2006;296:908–18. [36] Eftekhari SA. Dynamic analysis of a composite beam, rested on visco-elastic foundation, carrying multiple accelerating oscillators using a coupled finite element – differential quadrature method. MSc thesis, Shiraz University, Shiraz, 2006. [37] Eftekhari SA, Farid M. Dynamic analysis of multi-span laminated composite coated beams carrying multiple accelerating oscillators. ASME Appl Mech Dev 2006;257:83–92 [Paper No. IMECE2006-13872]. [38] Eftekhari SA, Farid M, Khani M. Dynamic analysis of laminated composite coated beams carrying multiple accelerating oscillators using a coupled finite element – differential quadrature method. ASME J Appl Mech 2009;76:1–13. [39] Eftekhari SA, Khani M. A coupled finite element–differential quadrature element method and its accuracy for moving load problem. Appl Math Model 2010;34:228–37. [40] Jafari AA, Eftekhari SA. A new mixed finite element–differential quadrature formulation for forced avibration of beams carrying moving loads. ASME J Appl Mech, submitted for publication. [41] Liu J, Wang X. An assessment of the differential quadrature time integration scheme for nonlinear dynamic equations. J Sound Vib 2008;314:246–53. [42] Civalek O. Harmonic differential quadrature-finite differences coupled approaches for geometrically nonlinear static and dynamic analysis of rectangular plates on elastic foundation. J Sound Vib 2006;294:966–80. [43] Civalek O. Nonlinear dynamic response of MDOF systems by the method of harmonic differential quadrature (HDQ). Int J Struct Eng Mech 2007;25(2): 201–17. [44] Civalek O. Nonlinear analysis of thin rectangular plates on Winkler–Pasternak elastic foundations by DSC–HDQ methods. Appl Math Model 2007;31:606–24. [45] Kim J, Paulino GH. Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. Int J Numer Methods Eng 2002;53:1903–35. [46] Zhang Ch, Savaidis A, Savaidis G, Zhu H. Transient dynamic analysis of a cracked functionally graded material by a BIEM. Comput Mater Sci 2003;26:167–74. [47] Wakashima K, Hirano T, Niino M. Space applications of advanced structural materials. ESP 1990;SP-303:397. [48] Reddy JN. Theory and analysis of elastic plates and shells. New York: CRC Press, Taylor & Francis Group; 2007. [49] Quan JR, Chang CT. New insights in solving distributed system equations by the quadrature methods. Part I. Analysis. Comput Chem Eng 1989;13:779–88. [50] Bathe KJ, Wilson EL. Numerical methods in finite element analysis. Englewood Cliffs (NJ): Prentice-Hall; 1976. [51] Meirovitch L. Analytical methods in vibrations. New York: Macmillan Company; 1967. [52] Rieker JR, Lin Y-H, Trethewey MW. Discretization considerations in moving load finite element beam models. Finite Element Anal Des 1996;21:129–44.