Dynamic behavior of nonuniform functionally graded euler-Bernoulli beams under multiple moving forces

Vietnam Journal of Mechanics, VAST, Vol. 37, No. 3 (2015), pp. 151 – 168 DOI:10.15625/0866-7136/37/3/4075 DYNAMIC BEHAVIOR OF NONUNIFORM FUNCTIONALLY GRADED EULER-BERNOULLI BEAMS UNDER MULTIPLE MOVING FORCES Le Thi Ha1, Nguyen Dinh Kien2,∗, Vu Tuan Anh3 1Hanoi University of Transport and Communications, Vietnam 2Institute of Mechanics, Vietnam Academy of Science and Technology, Hanoi, Vietnam 3Hanoi University of Science and Technology, Vietnam ∗E-mail: ndkien@imech.ac.vn Received June 02, 2014 Abstract. The dynamic behavior of nonuniform Euler-Bernoulli beams made of trans- versely functionally graded material under multiple moving forces is studied by the finite element method. The beam cross-section is assumed to vary in the width direction by two different types. A simple finite element formulation, accounting for variation of the mate- rial properties through the beam thickness and the shift in the physically neutral surface, is derived and employed in the study. The exact variation of the cross-sectional profile is employed in evaluation of the element stiffness and mass matrices. The dynamic response of the beam is computed with the aid of the implicit Newmark method. The numerical results show that the derived finite element formulation is capable to assess accurately the dynamic characteristics of the beam by using just several elements. The effect of the moving speed, material inhomogeneity and section profile on the dynamic behavior of the beams is investigated. The influence of the distance between the forces as well as the number of forces on the dynamic response is also examined and highlighted. Keywords: Functionally graded beam, physically neutral surface, moving force, dynamic behavior, finite element method. 1. INTRODUCTION Functionally graded materials (FGMs) have received much attention from engi- neers and researchers since they were first initiated by Japanese scientists in 1984 [1]. FGMs are produced by continuously varying volume fraction of constituent materials, usually ceramics and metals, in one or more spatial directions. As a result, the effective properties of FGMs exhibit continuous change, thus eliminating interface problems and mitigating thermal stress concentrations. Many investigations on analysis of FGM struc- tures subjected to different loadings are summarized in [2, 3], only contributions that are most relevant to the present work are briefly discussed below. c© 2015 Vietnam Academy of Science and Technology 152 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh Chakraborty et al. [4] proposed a first-order shear deformable beam element for analyzing the thermo-elastic behavior of FGM beams. In [5], the wave propagation be- havior of FGM beams under high frequency impulse loading was studied by using the spectral finite element method. Benatta et al. [6] derived an analytical solution to the bending problem of an FGM beam taking the warping effect into consideration. Based on the third-order shear deformation beam theory, Kadoli et al. [7] proposed a beam el- ement for studying the static behavior of FGM beams under ambient temperature. Lee et al. [8] presented a finite element procedure for computing the post-buckling response of FGM plates under compressive and thermal loads. Alshorbagy et al. [9], Shahba et al. [10, 11] derived beam finite elements for studying the free vibration of beams made of transversely and axially FGMs. Based on the concept of isogeometric analysis pro- posed by Hughes et al. [12], Tran et al. [13], Nguyen-Xuan et al. [14] developed the iso- geometric finite element formulations for static, dynamic and buckling analysis of FGM plates. The formulations utilized B-splines or non-uniform rational B-splines (NURBS) functions which enables to achieve easily the smoothness with arbitrary continuity or- der. Nguyen [15, 16], Nguyen and Gan [17] formulated nonlinear beam finite elements for investigating the large displacement behavior of tapered beams composed of axially and transversely FGMs. In [18] Nguyen et al. presented a finite element procedure for geometrically nonlinear analysis of planar FGM beam and frame structures. The problems of moving loads on an elastic beam are often met in the design of bridges, railways, highways. . . and they are subject of investigation for a long time. Both analytical method [16–18], and finite element method [19–22] are extensively employed in solving the moving load problems. With the rapid development and application of FGMs, analysis of FGM beams subjected to moving loads has been drawn attention from researchers recently. S¸ims¸ek and Kocatu¨rk [23] employed polynomials to approximate the displacement variables in solving the equations of motion of a transversely FGM Euler-Bernoulli beam subjected to a moving harmonic force. Also using the method in [23], S¸ims¸ek extended his work to problems of FGM beams subjected to a moving mass [24], and a nonlinear FGM beam under a moving harmonic force [25]. Rajabi et al. [26] studied the dynamic behavior of an FGM Euler-Bernoulli beam subjected to a moving oscillator by using the Runge-Kutta in solving the equations of motion. In [27], Malekzadeh and Monajjemzadeh used the finite element method to investigate the dy- namic response of an FGM plate resting on a Pasternak foundation subjected to thermal loading and a moving load. Also using the finite element method, Nguyen et al. [28] studied the vibration of a nonuniform FGM Timoshenko beam under a moving harmonic load. It has been stressed recently that the shift of the physically neutral surface of trans- versely FGM beams should be taken into account for correctly predicting the behavior of the beams [29]. In this line of work, Kang and Li [30, 31] determined the neutral axis position of a nonlinear FGM Euler-Bernoulli beam and then derived the solutions for tip displacements of the beam subjected to a tip moment or a tip transverse load. Based on the neutral surface and the third-order shear deformation beam theory, Zhang [32] in- vestigated the nonlinear bending of FGM beams. Eltaher et al. [33] considered the shift in the neutral axis position in the derivation of a beam finite element for studying the Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 153 free vibration of FGM macro/nano beams. It has been shown in [33] that the natural frequencies of an FGM beam are overestimated by ignoring the shift in the neutral axis position. Because of the interaction between the moving forces, the dynamic response of a beam to multiple moving forces, as shown by Henchi et al. [34], is very different from that of a beam subjected to a single moving force. To the authors’ best knowledge, the dynamic behavior of nonuniform FGM beams under multiple moving forces has not been studied in the literature, and it will be a subject of investigation of the present work. To this end, the finite element method previously used by the first two authors and their co- worker in Ref. [28], is again employed herein. The beam cross-section is assumed to vary in the width direction in two different manners. A finite element formulation, taking the variation of the elastic properties through the thickness and the shift in the physically neutral surface into account, is derived and employed in the study. It should be noted that in regard of the work by S¸ims¸ek and Kocatu¨rk [23], two different features are con- sidered in the present work. Firstly, the longitudinal variation of the beam cross-section is the one which is not easy to handle by the analysis method used in [23]. Secondly, the multiple moving forces, which has not been considered in [23] and in our previous work [28], requires some effort in numerical treatment. The dynamic response of the beam such as the time histories for mid-span deflection, dynamic deflection factor and axial stress distribution through the thickness are computed with the aid of the direct in- tegration Newmark method. The effect of the material inhomogeneity, section parameter and moving speed on the dynamic behavior of the beam is investigated in detail. The in- fluence of the material inhomogeneity, section profile and well as the loading parameters on the dynamic behavior of the beams is also examined and highlighted. 2. PROBLEM STATEMENT Fig. 1 shows a simply supported beam with length L, width b, height h, subjected to N forces P1, P2, ... PN , moving at a constant speed v from left to right. In the figure, a Cartesian co-ordinate system (x1, z1) is introduced as that the x1-axis lies on the bottom surface, and z1-axis directs upward. The distance between the force, d, is considered to be constant. The area A, and moment of inertia I of the beam cross-section are assumed to vary longitudinally in two following types - Type A: A = A0 ( 1− α ∣∣∣∣ xL − 12 ∣∣∣∣) , I = I0 (1− α ∣∣∣∣ xL − 12 ∣∣∣∣), - Type B: A = A0 [ 1− α ( x L − 1 2 )2] , I = I0 [ 1− α ( x L − 1 2 )2] , where A0 and I0 denote the area and moment of inertia of the mid-span cross-section, respectively; 0 ≤ α < 2 is the nonuniform section parameter. When α = 0, the beam becomes uniform. The two types of the section profile are depicted in the lower part of Fig. 1. The beam material is assumed to be composed of metal and ceramic phases whose volume fraction varies in the transverse direction according to 154 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh Vc = ( z1 h )n , Vc +Vm = 1, (1) where Vc, Vm are the volume fractions of ceramic and metal, respectively; n is the grading index, governing variation of the material properties through the beam thickness. As seen from Eq. (1), the bottom surface contains only metal and the top surface is pure ceramic. The composition is metal rich when n 1. Fig. 1. Nonuniform FGM beam under multiple moving forces 0 2 4 6 8 10 0.5 0.55 0.6 0.65 0.7 n h 0 /h E c /E m =1.5 E c /E m =2 E c /E m =3 E c /E m =5 Fig. 2. Dependence of neutral surface position on the index n The effective material properties (such as Young’s modulus and mass density), P , can be evaluated by a simple rule P(z1) = PcVc + PmVm = (Pc −Pm) ( z1 h )n + Pm, (2) where Pm and Pc are the properties of metal and ceramic, respectively. Clearly, due to the variation of the Young’s modulus E in the thickness direction, the neutral surface of the FGM beam does not coincide with the mid-plane. Denoting h0 as distance from the neutral surface to the bottom surface, and by introducing a new co-ordinate system (x, z) with the x-axis lies on the neutral surface, and z-axis directs upward as depicted in Fig. 1, the position of the neutral surface can be determined by using equilibrium condition for the beam subjected to pure bending as [29]∫ A σdA = b(x) ρ ∫ h−h0 −h0 zE(z)dz = 0, (3) where σ is the axial stress for the beam in pure bending, and ρ is the curvature radius of the neutral surface. Substituting z = (z1 − h0) into Eq. (3), the position of the neutral Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 155 surface can be determined with the aid of Eq. (2) as h0 = ∫ h 0 E(z1)z1dz1∫ h 0 E(z1)dz1 = h(n+ 1)(2Ec + nEm) 2(n+ 2)(Ec + nEm) . (4) The dependence of the neutral surface position upon the index n according to Eq. (4) is shown in Fig. 2 for various ratios of Young’s modulus of ceramic to that of metal, Ec/Em. As seen from the figure, for Ec/Em > 1 the physical neutral surface shifts upward from the mid-plane, regardless of the index n. Based on the Euler-Bernoulli beam theory, the displacements u1, u2, u3 at a point (x, y, z) in the x, y, z directions, respectively are given by u1(x, y, z, t) = u(x, t)− zw,x(x, t), u2(x, y, z, t) = 0, u3(x, y, z, t) = w(x, t), (5) where u(x, t) and w(x, t) are the axial and transverse displacements of a point on the neutral axis; z is a spatial co-ordinate in the thickness direction, and (...),x denotes the derivative with respect to x. Based on the Hook’s law, the axial strain e, and axial stress σ resulted from Eq. (5) are as follows e = u,x − zw,xx = u,x + zκ, σ = E(z)e = E(z)(u,x + zκ), (6) where κ = −w,xx is the beam curvature. The partial differential equations of motion for the beam under the moving forces can be derived by applying Hamilton’s principle. For the sake of brevity, the damping effect of the beam is not considered in the present work. The strain energy stored in the beam resulted from Eq. (6) has the following simple form U = 1 2 ∫ L 0 ∫ A(x) σe dAdx = 1 2 ∫ L 0 [ A11(x)u2,x − 2A12(x)u,xw,xx + A22(x)w2,xx ] dx, (7) in which (A11, A12, A22) = ∫ A(x) E(z)(1, z, z2)dA, (8) are the axial, axial-bending coupling and bending rigidities, respectively. It is worthy to note that in substituting z = z1 − h0 into Eq. (8), and taking Eq. (4) into consideration, the coupling rigidity A12 defined by Eq. (8) vanishes. As a result, the stiffness matrix resulted from Eq. (7) contains no coupling term. The kinetic energy of the beam resulted from the displacements (5) is as follows T = 1 2 ∫ L 0 ∫ A(x) ρ(z) ( u˙21 + u˙ 2 3 ) dAdx (9) = 1 2 ∫ L 0 [ I11(x)(u˙2 + w˙2)− 2I12(x)u˙w˙,x + I22(x)w˙2,x ] dx, 156 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh where a over dot indicates the derivative with respect to time t, and I11, I12 and I22 are the mass moments, defined as (I11, I12, I22) = ∫ A(x) ρ(z)(1, z, z2)dA, (10) where the mass density ρ(z) varies in the thickness direction according to Eq. (2) (with z = z1 − h0). It should be noted that for the longitudinal variation of cross-section con- sidered herein the rigidities Aij and the mass moments Iij depend upon x. In addition, the coupling mass moment I12, unfortunately does not vanish, and thus the mass matrix contains the coupling term. The potential energy of the moving forces is simply given by V = − N ∑ i=1 Piw(x, t)δ(x− vti(t)), (11) where δ(.) is the delta Dirac function, and ti is the time since the load Pi enters the beam from its left end. Applying Hamilton’s principle to Eqs. (7), (9) and (11), the differential equations of motion for the beam can be written in the forms I11u¨− I12w¨,x − (A11u,x),x = 0, I11w¨+ (I12u¨),x − (I22w¨,x),x + (A22w,xx),xx = N ∑ i=1 Piδ(x− vti). (12) Except for the presence of the coupling mass moment I12, the system of equations (12) has the same forms as that of a nouniform homogeneous beam. 3. FINITE ELEMENT FORMULATION The finite element method is employed herein to solve Eq. (12). To this end, the beam is assumed being divided into a number of two-node beam elements with length of l. There are axial and transverse displacements and a rotation at each node. Thus, the vector of nodal displacements, d, for a generic element has the following components d = {ui wi θi uj wj θj}T, (13) where and hereafter a superscript ‘T’ denotes the transpose of a vector or a matrix. The axial displacement u and transverse displacement w are interpolated from the nodal dis- placements according to u = Nud , w = Nwd, (14) where Nu and Nw are the matrices of shape functions for u and w, respectively. Substi- tuting Eq. (14) into Eqs. (7) and (9), we get U = 1 2 nel ∑ i=1 dTkd = 1 2 nel ∑ i=1 dT(kuu + kθθ)d, (15) and T = 1 2 nel ∑ i=1 d˙Tmd˙ = 1 2 nel ∑ i=1 d˙T(muu + mww + muθ + mθθ)d˙. (16) Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 157 In Eqs. (15) and (16), nel is the total number of elements; k and m are respectively the element stiffness and mass matrices, and kuu = ∫ l 0 NTu,xA11Nu,xdx, kθθ = ∫ l 0 NTw,xxA22Nw,xxdx, (17) are respectively the stiffness matrices stemming from stretching and bending, muu = ∫ l 0 NTu I11Nudx , mww = ∫ l 0 NTw I11Nwdx muθ = ∫ l 0 NTu I12Nw,xdx , mθθ = ∫ l 0 NTw,x I22Nw,xdx (18) are the mass matrices stemming from axial displacement, transverse displacement, axial- bending coupling and cross-section rotation, respectively. Having the element stiffness and mass matrices derived, the finite element equa- tion for vibration of the beam is as follows MD¨ + KD = Fex, (19) where M, K are the structural mass and stiffness matrices assembled from the element mass and stiffness matrices, respectively; Fex is the structural nodal load vector of the external forces with the following form Fex = 0 ... P1Nw|x1︸ ︷︷ ︸ loading element 0 ... 0 PiNw|xi︸ ︷︷ ︸ loading element 0 ... PNNw|xN︸ ︷︷ ︸ loading element 0 ... 0  T , (20) which contains all zero coefficients, except for the elements currently under loading. The notation NTw|xi in the above equation implies that the shape functions Nw are evaluated at the abscissa xi, the current position of load Pi. The system of equations (19) can be solved by the direct integration Newmark method. The average acceleration implicit Newmark method described in [35], which ensures the unconditional convergency is adopted in the present work. In the free vi- bration analysis, the right hand side of Eq. (19) is set to zeros, and a harmonic response, D = D¯ sinωt is assumed, so that Eq. (19) deduces to (K−ω2M)D¯ = 0, (21) where ω is the circular frequency, and D¯ is the vibration amplitude. Eq. (21) can be solved by a standard method of the eigenvalue problem [35]. To improve the accuracy, the exact variation of the cross-sectional profiles is employed in evaluation of the rigidities and mass moments defined in Eqs. (8) and (10), respectively. 158 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh 4. NUMERICAL RESULTS AND DISCUSSION A simply supported FGM beam with L = 20 m, h = 0.8 m, b0 = 2 m, where b0 is the width of the mid-span cross section, is employed in this Section to study the dynamic response of the beam. Steel and alumina are employed as metal and ceramic phases of the FGM, respectively. The Young’s modulus and mass density are respectively Em = 210 GPa and ρm = 7800 kg/m3 for steel, and that for alumina are Ec = 390 GPa and ρc = 3960 kg/m3 [23]. Unless stated, the beam is assumed under action of three moving forces with the same amplitude, P1 = P2 = P3 = P0 = 100 kN. Linear and cubic Hermite polynomials are adopted as the shape functions for the axial and transverse displacements, respectively. Thus, the matrices of shape functions Nu and Nw in Eq. (14) have the following forms Nu = {Nu1 0 0 Nu2 0 0}, Nw = {0 Nw1 Nw2 0 Nw3 Nw4}, (22) in which Nu1 = l − x l , Nu2 = x l , Nw1 = 2 x3 l3 − 3x 2 l2 + 1 , Nw2 = x3 l2 − 2x 2 l + x , Nw3 = −2x 3 l3 + 3 x2 l2 , Nw4 = x3 l2 − 2x 2 l . (23) For the case of constant moving speed considered herein, total time ∆T necessary for a force to cross the beam is L/v. In the computation reported below a uniform time incre- ment width of ∆t = ∆T/500 is employed for the Newmark method. In order to facilitate the discussion of numerical results, the following dimensionless parameters represent- ing the maximum mid-span dynamic deflection and the moving force speed are intro- duced as fD = max ( w(L/2, t) w0 ) , fv = v v0cr = piv ω01L , (24) where w0 = P0L3/48Em I0 is the static deflection of the uniform steel beam under a static load P0 acting at the mid-span; v0cr = ω01L/pi, with ω 0 1 = pi2 L2 √ Em I0/ρmA0, is the critical speed of the simply supported uniform steel beam [19]. The definition of parameter fD by Eq. (24) is similar to that of the dynamic magnification factor in the moving load problem of homogeneous beams [19]. However, for the FGM beam considered in the present work, fD is not only governed by the moving speed but by the material inhomogeneity and the section profile also, and it will be called the dynamic deflection factor in the below. 4.1. Formulation verification In order to verify the accuracy of the derived finite element formulation, the fun- damental frequency and dynamic response of a uniform FGM Bernoulli beam subjected to a moving point force are firstly computed and compared to the result of Ref. [23]. To this end, a simply supported beam with width b = 0.4 m, height h = 0.9 m, previously Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 159 used in Ref. [23], is adopted in this subsection. The beam is also composed of steel and alumina with the above mentioned material properties. In Tab. 1, the fundamental frequency parameter of a uniform FGM beam with as- sumed properties Ec/Em = 3, ρc = ρm obtained by different numbers of elements is listed for various values of the index n and different length to height ratios, L/h = 20 and L/h = 100. The frequency parameter µ1 in Tab. 1 is defined as µ1 = ω1L2 √ ρmA/Em I, where A = bh, I = bh3/12, and ω1 is the fundamental frequency of the beam. The Tab. 1 shows the fast convergency of the present formulation, and all the frequencies converge by using just ten elements. In Tab. 2, the fundamental frequency parameter is given for various values of the index n, the ratio of length to height L/h, and the ratio of Young’s moduli Erat = Ec/Em. Due to the convergency stated above, only twelve elements were used in evaluating the frequencies in Tab. 2. The corresponding parameter obtained by S¸ims¸ek and Kocatu¨rk in Ref. [23] is also listed in Tab. 2. Tab. 2 shows the good agreement between the fundamental frequencies obtained in the present work with that of Ref. [23]. In Tab. 3, the maximum dynamic deflection factor and the corresponding moving speed of the uniform FGM beam are listed for various values of the index n. For comparison purpose, the corresponding data of Ref. [23] are also given in Tab. 3. Very good agree- ment between the numerical result of the present work with that of Ref. [23] is noted. Table 1. Convergency of present formulation in evaluating fundamental frequency of a uniform FGM beam (ρc = ρm, Ec/Em = 3) nel L/h n 2 4 6 8 10 12 Ref. [23] 20 0.1 4.0555 4.0481 4.0477 4.0476 4.0476 4.0476 4.0475 0.2 3.9820 3.9747 3.9742 3.9742 3.9741 3.9741 3.9741 2 3.5386 3.5321 3.5317 3.5317 3.5317 3.5317 3.5308 3 3.4935 3.4871 3.4867 3.4867 3.4867 3.4867 3.4858 10 3.3810 3.3748 3.3745 3.3744 3.3744 3.3744 3.3738 100 0.1 4.0572 4.0497 4.0493 4.0492 4.0492 4.0492 4.0495 0.2 3.9836 3.9763 3.9758 3.9758 3.9758 3.9758 3.9761 2 3.5402 3.5337 3.5333 3.5333 3.5333 3.5333 3.5331 3 3.4951 3.4887 3.4883 3.4882 3.4882 3.4882 3.4881 10 3.3825 3.3762 3.3759 3.3758 3.3758 3.3758 3.3757 Secondary, the time history for dynamic mid-span deflection of a uniform homo- geneous beam subjected to three moving forces, previously studied by Henchi et al. in Ref. [34] by the dynamic stiffness method, is computed. The beam geometric and ma- terial data are: L = 24.384 m, A = 0.954 m2, I = 2.9× 10−4 m4, E = 19× 1011 N/m2, 160 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh Table 2. Fundamental frequency parameter µ1 of uniform FGM beam with assumed material properties ρc = ρm (Erat = Ea/Es) L/h Erat Source n = 0.1 n = 0.2 n = 1 n = 2 n = 3 n = 10 20 2 Ref. [23] 3.6775 3.6301 3.4421 3.3765 3.3500 3.2725 Present 3.6776 3.6303 3.4426 3.3770 3.3505 3.2729 4 Ref. [23] 4.3370 4.2459 3.8234 3.6485 3.5858 3.4543 Present 4.3370 4.2459 3.8243 3.6496 3.5870 3.4551 100 2 Ref. [23] 3.6793 3.6320 3.4440 3.3784 3.3519 3.2742 Present 3.6791 3.6318 3.4440 3.3784 3.3519 3.2743 4 Ref. [23] 4.3392 4.2481 3.8259 3.6513 3.5886 3.4565 Present 4.3388 4.2476 3.8260 3.6514 3.5887 3.4566 Table 3. Maximum dynamic deflection factor and corresponding moving speed of uniform FGM beam under a single moving force n max( fD) v (m/s) Present Ref. [23] Present Ref. [23] 0.2 1.0347 1.0344 222 222 0.5 1.1445 1.1444 197 198 1 1.2504 1.2503 179 179 2 1.3377 1.3376 164 164 pure steel 1.7326 1.7324 132 132 pure alumina 0.9329 0.9328 252 252 0 0.4 0.8 1.2 1.6 2 −0.5 0 0.5 1 1.5 2 2.5 3 t/ΔT w (L/ 2,t ) ( mm /0. 30 48 ) Henchi et al. (1997) present work Fig. 3. Dynamic mid-span deflection of uniform homogeneous beam under three moving forces (d = L/4, v = 22.5 m/s) Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 161 ρA = 9.576 × 103 kg/m3, where L, A, I, E, ρ are the total length, cross-sectional area, moment of inertia of cross-section, Young’s modulus and mass density of the beam, re- spectively. Fig. 3 shows the mid-span dynamic deflection for the case P0 = 5324.256 N, v = 22.5 m/s, d = L/4, where the numerical result obtained by Henchi et al. [34] is also depicted. The figure shows a good agreement between the finite element solution of the present work with the result obtained by the dynamic stiffness matrix method of Ref. [34]. 4.2. Effect of material inhomogeneity In Fig. 4 the time histories for mid-span deflection of the type A beam with α = 0.5 are depicted for various values of the index n and two values of the speed parameter, fv = 1/8 and fv = 1/4. At a given value of the moving speed, the dynamic deflection of the beam, as seen from the figure, is greatly influenced by the material parameter n. The maximum dynamic deflection of the beam steadily increases when rasing the index n, regardless of the moving speed. This can be explained by the fact that, as seen from Eq. (1), the beam with a higher index n contains more steel, and thus it is softer. The increase in the maximum mid-span dynamic deflection is also clearly seen from Fig. 5, where the deflection factor fD is shown as a function of the speed parameter fv for various values of the index n. 0 0.5 1 1.5 −0.5 0 0.5 1 1.5 2 2.5 t/ΔT w (L/ 2,t )/w 0 0 0.5 1 1.5 −0.5 0 0.5 1 1.5 2 2.5 t/ΔT w (L/ 2,t )/w 0 n=0.2 n=0.5 n=2 n=5 n=0.2 n=0.5 n=2 n=5(a) f v =1/8 (b) fv=1/4 Fig. 4. Time histories for mid-span deflection of type A beam under three moving forces (α = 0.5, d = L/4) In Fig. 6, the distribution through the beam thickness of the axial stress at the mid- span section is depicted for different values of the index n and two values of the speed parameter, fv = 1/8 and fv = 1/4. The axial stress shown in the figure was normalized by the maximum static axial stress of a uniform steel beam, σ0 = P0Lh/8I0, and it was computed at the time when the second force arrives at the midpoint of the beam. As seen from the figure, the axial stress distribution of the FGM beam is very different from that of the homogeneous beam. The stress of the FGM beam is not symmetrical with regard to the coordinate original, and it does not become zero on the mid-plane, regardless of 162 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh 0 0.5 1 1.5 2 1.5 2 2.5 3 3.5 f v f D n=0.2 n=0.5 n=2 n=5 Fig. 5. Speed parameter versus deflection factor of type A beam subjected to three moving forces (α = 0.5, d = L/4) −4 −3 −2 −1 0 1 2 3 −0.5 −0.25 0 0.25 0.5 σ/σ0 z/ h −4 −3 −2 −1 0 1 2 3 −0.5 −0.25 0 0.25 0.5 σ/σ0 z/ h pure steel n=0.5 n=5 pure steel n=0.5 n=5 (a) f v =1/8 (b) f v =1/4 Fig. 6. Normalized axial stress distribution through the thickness of type A beam subjected to three moving forces (α = 0.5, d = L/8) the moving speed. The moving speed slightly alters the amplitude of the stress, but it hardly changes the distribution of the stress. 4.3. Effect of moving speed The influence of the moving speed on the dynamic deflection factor of the FGM beam is clearly seen from Fig. 5. For lower values of the speed parameter fv, the deflec- tion factor in Fig. 5 both increases and decreases with increasing fv, and this phenomenon is associated with the oscillations of the beam when it subjected to a low speed moving load [36]. For higher values of fv, as in case of homogeneous beams, the deflection factor increases when raising the speed parameter fv, it then reaches a maximum value before decreases. The effect of the moving speed can also be seen from the time histories for the Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 163 0 0.5 1 1.5 −1 0 1 2 3 4 t/ΔT w (L/ 2,t )/w 0 0 0.5 1 1.5 −1 0 1 2 3 4 t/ΔT w (L/ 2,t )/w 0 f v =1/8 f v =1/4 f v =1/2 f v =2/3 f v =1/8 f v =1/4 f v =1/2 f v =2/3(a) α=0.5 (b) α=1.5 Fig. 7. Time histories for mid-span deflection of type A beam under different speeds moving point forces (n = 3, d = L/4) 0 2 4 6 8 10 1.5 2 2.5 3 3.5 4 n f D 0 2 4 6 8 10 1.5 2 2.5 3 3.5 4 n f D f v =1/8 f v =1/4 f v =1/2 f v =1/8 f v =1/4 f v =1/2 (a) α=0.5 (b) α=1.5 Fig. 8. Material index n versus deflection factor fD of type A beam subjected to three moving forces (d = L/4) mid-span deflection and the relation between the deflection factor fD and the index n as depicted in Fig. 7 and Fig. 8 for the type A beam, respectively. 4.4. Effect of distance between the forces In Fig. 9, the relation between the deflection factor fD and the moving speed pa- rameter fv of the type A beam with n = 0.5 is shown for various values of the distance between the forces d and the section parameter α. The effect of the distance between the forces is clearly seen from the figure, where the deflection factor fD is remarkably lower for a smaller distance d, regardless of the speed and section parameters. The similar sit- uation is occurred for the axial stress as depicted in Fig. 10, where the distribution of the stress through the beam thickness is shown for the type A beam with an index n = 0.5. 164 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh 0 0.5 1 1.5 2 0.5 1 1.5 2 2.5 3 3.5 4 f v f D 0 0.5 1 1.5 2 0.5 1 1.5 2 2.5 3 3.5 4 f v f D d=L/8 d=L/4 d=L/2 d=L/8 d=L/4 d=L/2(a) α=0.5 (b) α=1 Fig. 9. Effect of distance between moving forces on relation between deflection factor and speed parameter of type A beam (n = 0.5) −3 −2 −1 0 1 2 −0.5 −0.25 0 0.25 0.5 σ x /σ0 z/h 0 −3 −2 −1 0 1 2 −0.5 −0.25 0 0.25 0.5 σ x /σ0 z/h 0 d=L/8 d=L/4 d=L/2 d=L/8 d=L/4 d=L/2 (a) f v =1/8 (b) fv=1/4 Fig. 10. Effect of distance between moving forces on thickness distribution of axial stress of type A beam (n = 0.5) The axial stress increases considerably when the distance between the forces is smaller, regardless of the moving speed. 4.5. Effect of section profile In Tab. 4 the maximum values of the dynamic deflection factor of the FGM beam subjected to three moving forces are listed for the two types of the section profile, and for various values of the section parameter α and the index n. The maximum deflection factor listed in the table increases by raising the section parameter α, regardless of the section type and the index n. The maximum deflection factor of the type A beam, as seen from the table, is much more sensitive to the change in the section parameter α compares to that of the type B beam. For example, with n = 3 the maximum deflection factor of the type A beam increases 18.07% when raising the section parameter from 0 to 1.2, while this value is just 3.54% for the type B beam. The effect of the section profile on the dynamic response of the FGM beam can also be seen from the relation between the Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 165 deflection factor fD and the index n as shown in Fig. 11 for two values of the distance between the forces, d = L/8 and d = L/4. The sensitivity of the type A beam with the change in the section parameter is again clearly observed from the figure. Table 4. Maximum dynamic deflection factor, max( fD), for FGM beam with different section profiles subjected to three moving forces (d = L/8) Section n 0 0.2 0.4 0.6 0.8 1 1.2 Type A 0.2 2.8729 2.9535 3.0417 3.1388 3.2467 3.3680 3.5064 0.5 3.1776 3.2668 3.3643 3.4717 3.5911 3.7253 3.8784 3 3.8203 3.9276 4.0448 4.1740 4.3175 4.4788 4.6629 5 3.9509 4.0618 4.1830 4.3166 4.4650 4.6318 4.8222 Type B 0.2 2.8729 2.8891 2.9058 2.9231 2.9409 2.9594 2.9785 0.5 3.1956 3.2140 3.2331 3.2528 3.2732 3.2944 3.3163 3 3.8203 3.8419 3.8642 3.8871 3.9108 3.9353 3.9607 5 3.9509 3.9732 3.9962 4.0199 4.0444 4.0698 4.0961 0 2 4 6 8 10 1.5 2 2.5 3 3.5 n f D 0 2 4 6 8 10 1 1.5 2 2.5 3 3.5 n f D type A, α=0.5 type B, α=0.5 type A, α=1.5 type B, α=1.5 type A, α=0.5 type B, α=0.5 type A, α=1.5 type B, α=1.5 (a) d=L/8 (b) d=L/4 Fig. 11. Relation between deflection factor and material index of FGM beam with different section profiles under three moving forces ( fv = 0.5) 4.6. Effect of different force numbers The time histories for the mid-span deflection of type A beam subjected to different numbers of moving forces are shown in Fig. 12 for n = 3, α = 0.5, d = L/4, and for two values of the speed parameter, fv = 1/8 and fv = 1/4. As expected, the maximum mid-span dynamic deflection, as seen from Fig. 12, increases when the beam subjected to more numbers of the moving forces. The number of the forces also changes the time at which the maximum mid-span deflection occurs. The relation between the dynamic deflection factor and the moving speed parameter depicted in Fig. 13 for the type A beam with n = 3 and α = 0.5 under different numbers of the moving forces also clearly shows the increase in the deflection factor when the beam subjected to more numbers of the moving forces. 166 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh 0 0.5 1 1.5 2 −0.5 0 0.5 1 1.5 2 2.5 t/∆T w(L /2,t )/w 0 0 0.5 1 1.5 2 −0.5 0 0.5 1 1.5 2 2.5 t/∆T w(L /2,t )/w 0 (a) f v =1/8 (b) f v =1/4 1 2 3 4 N: 1 2 3 4 N:  ó k k  C ó k ó  Fig. 12. Time histories for mid-span deflection of type A beam subjected to different numbers of moving forces (n = 3, α = 0.5, d = L/4) ó 0 0.5 1 1.5 1 2 3 4 5 6 f v f D 0 0.5 1 1.5 1 2 3 4 f v f D (a ) d=L/8 (b ) d=L/4 N: 4 3 2 4 3 2 N:  C ó k ó  Fig. 13. Relation between deflection factor and speed parameter of type A beam subjected to different numbers of moving forces (n = 3, α = 0.5) 5. CONCLUSIONS The dynamic behavior of nonuniform FGM Euler-Bernoulli beams subjected to multiple moving forces has been studied by using the finite element method. The ma- terial properties of the beams are assumed to vary in the thickness direction by a power law function. A finite element formulation, taking the effect of the cross-section varia- tion and the material inhomogeneity into account, has been derived and employed in the study. The exact variation of the section profile was used in evaluation of the ele- ment formulation. The dynamic response of the beam was computed with the aid of the implicit Newmark method. The numerical results have shown that the derived element formulation is accurate in evaluating the dynamic response of the beams. The dynamic characteristics, including the time history, dynamic deflection factor, axial stress distribu- tion are governed by the moving speed, section profile, number of forces as well as the distance between the forces. A parametric study has been carried out to highlight the influence of the material inhomogeneity, section profile and loading parameters on the dynamic behavior of the nonuniform FGM beams under the movings forces. Dynamic behavior of nonuniform functionally graded Euler-Bernoulli beams under multiple moving forces 167 REFERENCES [1] M. Koizumi. FGM activities in Japan. Composites Part B: Engineering, 28, (1), (1997), pp. 1–4. [2] V. Birman and L. W. Byrd. Modeling and analysis of functionally graded materials and struc- tures. Applied Mechanics Reviews, 60, (5), (2007), pp. 195–216. [3] D. K. Jha, T. Kant, and R. K. Singh. A critical review of recent research on functionally graded plates. Composite Structures, 96, (2013), pp. 833–849. [4] A. Chakraborty, S. Gopalakrishnan, and J. N. Reddy. A new beam finite element for the analysis of functionally graded materials. International Journal of Mechanical Sciences, 45, (3), (2003), pp. 519–539. [5] A. Chakraborty and S. Gopalakrishnan. A spectrally formulated finite element for wave propagation analysis in functionally graded beams. International Journal of Solids and Struc- tures, 40, (10), (2003), pp. 2421–2448. [6] M. A. Benatta, I. Mechab, A. Tounsi, and E. A. A. Bedia. Static analysis of functionally graded short beams including warping and shear deformation effects. Computational Materials Sci- ence, 44, (2), (2008), pp. 765–773. [7] R. Kadoli, K. Akhtar, and N. Ganesan. Static analysis of functionally graded beams us- ing higher order shear deformation theory. Applied Mathematical Modelling, 32, (12), (2008), pp. 2509–2525. [8] Y. Y. Lee, X. Zhao, and J. N. Reddy. Postbuckling analysis of functionally graded plates sub- ject to compressive and thermal loads. Computer Methods in Applied Mechanics and Engineer- ing, 199, (25), (2010), pp. 1645–1653. [9] A. E. Alshorbagy, M. A. Eltaher, and F. F. Mahmoud. Free vibration characteristics of a functionally graded beam by finite element method. Applied Mathematical Modelling, 35, (1), (2011), pp. 412–425. [10] A. Shahba, R. Attarnejad, and S. Hajilar. Free vibration and stability of axially functionally graded tapered Euler-Bernoulli beams. Shock and Vibration, 18, (5), (2011), pp. 683–696. [11] A. Shahba, R. Attarnejad, M. T. Marvi, and S. Hajilar. Free vibration and stability analysis of axially functionally graded tapered timoshenko beams with classical and non-classical boundary conditions. Composites Part B: Engineering, 42, (4), (2011), pp. 801–808. [12] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194, (39), (2005), pp. 4135–4195. [13] L. V. Tran, A. J. M. Ferreira, and H. Nguyen-Xuan. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites Part B: Engineering, 51, (2013), pp. 368–383. [14] H. Nguyen-Xuan, L. V. Tran, C. H. Thai, S. Kulasegaram, and S. P. A. Bordas. Isogeometric analysis of functionally graded plates using a refined plate theory. Composites Part B: Engi- neering, 64, (2014), pp. 222–234. [15] D. K. Nguyen. Large displacement response of tapered cantilever beams made of axially functionally graded material. Composites Part B: Engineering, 55, (2013), pp. 298–305. [16] D. K. Nguyen. Large displacement behaviour of tapered cantilever Euler-Bernoulli beams made of functionally graded material. Applied Mathematics and Computation, 237, (2014), pp. 340–355. [17] D. K. Nguyen and B. S. Gan. Large deflections of tapered functionally graded beams sub- jected to end forces. Applied Mathematical Modelling, 38, (11), (2014), pp. 3054–3066. 168 Le Thi Ha, Nguyen Dinh Kien, Vu Tuan Anh [18] D. K. Nguyen, B. S. Gan, and T. H. Trinh. Geometrically nonlinear analysis of planar beam and frame structures made of functionally graded material. Structural Engineering and Me- chanics, 49, (6), (2014), pp. 727–743. [19] L. Fry`ba. Vibration of solids and structures under moving loads. Academia, Prague, (1972). [20] S. Chonan. The elastically supported Timoshenko beam subjected to an axial force and a moving load. International Journal of Mechanical Sciences, 17, (9), (1975), pp. 573–581. [21] H. P. Lee. Dynamic response of a beam with intermediate point constraints subject to a mov- ing load. Journal of Sound and Vibration, 171, (3), (1994), pp. 361–368. [22] Y. H. Lin and M. W. Trethewey. Finite element analysis of elastic beams subjected to moving dynamic loads. Journal of Sound and Vibration, 136, (2), (1990), pp. 323–342. [23] M. S¸ims¸ek and T. Kocatu¨rk. Free and forced vibration of a functionally graded beam sub- jected to a concentrated moving harmonic load. Composite Structures, 90, (4), (2009), pp. 465– 473. [24] M. S¸ims¸ek. Vibration analysis of a functionally graded beam under a moving mass by using different beam theories. Composite Structures, 92, (4), (2010), pp. 904–917. [25] M. S¸ims¸ek. Non-linear vibration analysis of a functionally graded Timoshenko beam under action of a moving harmonic load. Composite Structures, 92, (10), (2010), pp. 2532–2546. [26] K. Rajabi, M. H. Kargarnovin, and M. Gharini. Dynamic analysis of a functionally graded simply supported Euler-Bernoulli beam subjected to a moving oscillator. Acta Mechanica, 224, (2), (2013), pp. 425–446. [27] P. Malekzadeh and S. M. Monajjemzadeh. Dynamic response of functionally graded plates in thermal environment under moving load. Composites Part B: Engineering, 45, (1), (2013), pp. 1521–1533. [28] D. K. Nguyen, B. S. Gan, and T. H. Le. Dynamic response of non-uniform functionally graded beams subjected to a variable speed moving load. Journal of Computational Science and Tech- nology, 7, (1), (2013), pp. 12–27. [29] S. V. Levyakov. Elastica solution for thermal bending of a functionally graded beam. Acta Mechanica, 224, (8), (2013), pp. 1731–1740. [30] Y.-A. Kang and X.-F. Li. Bending of functionally graded cantilever beam with power-law non-linearity subjected to an end force. International Journal of Non-Linear Mechanics, 44, (6), (2009), pp. 696–703. [31] Y.-A. Kang and X.-F. Li. Large deflections of a non-linear cantilever functionally graded beam. Journal of Reinforced Plastics and Composites, 29, (12), (2010), pp. 1761–1774. [32] D.-G. Zhang. Nonlinear bending analysis of FGM beams based on physical neutral surface and high order shear deformation theory. Composite Structures, 100, (2013), pp. 121–126. [33] M. A. Eltaher, A. E. Alshorbagy, and F. F. Mahmoud. Determination of neutral axis position and its effect on natural frequencies of functionally graded macro/nanobeams. Composite Structures, 99, (2013), pp. 193–201. [34] K. Henchi, M. Fafard, G. Dhatt, and M. Talbot. Dynamic behaviour of multi-span beams under moving loads. Journal of Sound and Vibration, 199, (1), (1997), pp. 33–50. [35] M. Ge´radin and D. J. Rixen. Mechanical vibrations: Theory and application to structural dynamics. John Wiley & Sons, 2nd edition, (1997). [36] M. Olsson. On the fundamental moving load problem. Journal of Sound and Vibration, 145, (2), (1991), pp. 299–307.