Dynamic Behavior of a Bidirectional Functionally Graded Sandwich Beam under Nonuniform Motion of a Moving Load

Abstract

A bidirectional functionally graded Sandwich (BFGSW) beam model made from three distinct materials is proposed and its dynamic behavior due to nonuniform motion of a moving point load is investigated for the first time. The beam consists of three layers, a homogeneous core, and two functionally graded face sheets with material properties varying in both the thickness and longitudinal directions by power gradation laws. Based on the first-order shear deformation beam theory, a finite beam element is derived and employed in computing dynamic response of the beam. The element which used the shear correction factor is simple with the stiffness and mass matrices evaluated analytically. The numerical result reveals that the material distribution plays an important role in the dynamic response of the beam, and the beam can be designed to meet the desired dynamic magnification factor by appropriately choosing the material grading indexes. A parametric study is carried out to highlight the effects of the material distribution, the beam layer thickness and aspect ratios, and the moving load speed on the dynamic characteristics. The influence of acceleration and deceleration of the moving load on the dynamic behavior of the beam is also examined and highlighted.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Investigations on vibration of beams made of functionally graded materials (FGMs), a new type of composite materials initiated by Japanese researchers in 1984 [1], have been extensively carried out in the last two decades. It has been shown that the natural frequencies of FGM beams are governed by the variation of the material properties in the beam thickness or length direction [2–6]. Dynamic behavior of FGM beams under moving loads, the topic discussed herein, has also been reported in recent years. In this line of works, Şimşek and Kocatürk [7] and Şimşek [8] employed polynomials to approximate the displacement field in computing the dynamic response of FGM beams to moving loads. The authors showed that the dynamic deflection of the beams is significantly influenced by the material gradation in the beam thickness. The dynamic behavior of beams due to a moving harmonic load is also significantly influenced by the variation of the material properties in the beam length [9]. The Rayleigh-Ritz method was used in combination with the differential quadrature method by Khalili et al. [10] in examining dynamic behavior of an Euler-Bernoulli beam due to a moving mass. The material properties are considered to follow an exponential or a power-law function. The forced vibration of an Euler-Bernoulli beam under a moving oscillator was studied by Rajabi et al. [11] using the Petrov-Galerkin method. Gan et al. [12] formulated a finite element formulation for computing dynamic response of a Timoshenko beam with material properties varying along the beam length. Linear variation of the beam width was considered by the authors. The dynamic response of an axially FGM beam to a moving harmonic force was examined by Wang and Wu [13] using the Lagrange method. The effect of axial compressive force and the temperature rise were taken into consideration by the authors. Nguyen and Bui [14] studied the effect of temperature rise on dynamic response of a FGM Timoshenko beam under a moving load using a hierarchical beam element. Finite element method was also used by Esen et al. [15] in computing the dynamic response of a FGM Timoshenko to an accelerating moving mass.

With the development of the advanced manufacturing methods [16], FGMs can be employed to fabricate structural Sandwich elements to improve performance of structures. FGM Sandwich structures with smooth variation of material properties can eliminate the interface separation problems which are often seen in the conventional Sandwich structures. Analysis of FGM Sandwich beams has drawn much attention from researchers recently. Chakraborty et al. [17] presented a finite element procedure for thermoelastic analysis of FGM and FGM Sandwich Timoshenko beams. Based on different beam models, Apetre et al. [18] examined the bending behavior of Sandwich beams with a FGM core. Adopting a higher-order sandwich panel theory, Rahmani et al. [19] studied free vibration of Sandwich beams with a syntactic core. Pradhan and Murmu [20] used the modified differential quadrature method to study thermomechanical vibration of FGM Sandwich beams resting on an elastic foundation. A refined sinusoidal shear deformation beam theory was used by Zenkour et al. [21] to study static bending of a FGM Sandwich beam resting on a Pasternak foundation. Su et al. [22] employed the general Fourier formulation to calculate frequencies of FGM Sandwich beams on an elastic foundation. Free vibration analysis of FGM Sandwich beams was carried out by Amirani et al. [23] and Yang et al. [24] using the element-free Galerkin and mesh-free radial point interpolation methods. Vo et al. [25–27] and Nguyen et al. [28–30] presented various higher-order shear deformation theories for buckling, free vibration, and bending analyses of FGSW beams. The transverse displacement in the theories is split into bending and shear parts, and the effect of the thickness stretching is taken into account. The bucking and free vibration of FGM Sandwich beams were also considered in [31] using a refined hyperbolic shear and normal deformation beam theory. Lagrange multiplier method was used in combination with Newmark method by Şimşek and Al-shujairi [32] to evaluate dynamic response of FGM Sandwich beams under two moving forces. Songsuwan et al. [33] employed the Ritz and Newmark methods to study dynamic behavior of FGM Sandwich beams under a moving harmonic load.

Development of structural elements with material properties varying in two or more directions for withstanding severe general loadings is of great important in practice. Several bidirectional FGM beam models and their mechanical behavior have been considered in recent years. Şimşek [34] and Hao and Wei [35] assumed an exponential variation for the material properties in both the thickness and length directions in their vibration analysis of FGM beams. The authors showed that the dynamic characteristics of the beams are greatly influenced by the material gradient distribution, and the material properties of the bidirectional FGM beams can be designed to meet the goals of optimization. Wang et al. [36] considered the material properties varying by an exponential function along the beam length and a power law though the thickness in their free vibration study of bidirectional FGM Euler-Bernoulli beams. Free vibration analysis of FGM beams with material properties varying in the length and thickness directions by various gradation laws was carried out by Huynh et al. [37] using the NURBS isogeometric finite element approach. Based on a quasi-3D beam theory, Karamanli [38] studied static bending of a FGM Sandwich beam with material properties varying in the length and thickness direction by the power laws. The symmetric smoothed particle hydrodynamics (SSPH) method was employed to evaluate the response of the beam to a distributed load. Finite element method was used by Nguyen et al. [39] and Nguyen and Tran [40] to study vibration of bidirectional FGM Timoshenko beams. The beams were considered to be made from four materials with volume fraction varying in both the thickness and length directions by the power gradation laws.

It is clear from the above literature review that, except for the work in [38], the mechanical behavior of bidirectional functionally graded Sandwich (BFGSW) beams has not been studied so far. This paper tries to make an effort to address this shortage. To this end, a BFGSW beam model is proposed and its dynamic behavior due to nonuniform motion of a moving point load is studied. The beam made from three distinct materials has three layers, a homogeneous core, and two FGM face sheets with material properties varying in both the thickness and length directions by power gradation laws. The conventional unidirectional FGM Sandwich beam, for example, the Sandwich beam in [32, 33], is a special case of the present beam model. In addition to the BFGSW beam model, the main novelty of this paper is that the effect of variation of the material properties in both the thickness and longitudinal directions on the dynamic behavior of FGM Sandwich beams under moving loads is taken into consideration herein for the first time. Based on the first-order shear deformation beam theory, a finite beam element is formulated and employed in computing the dynamic response of the beam. The first-order shear deformation beam theory requires a shear correction factor [41, 42], but it leads to a simple element. The element with stiffness and mass matrices evaluated analytically is derived using quadratic and cubic polynomials to interpolate the rotation and transverse displacement, respectively. The effects of the material distribution, the layer thickness and aspect ratios, and the moving load parameters on the dynamic behavior of the beam are investigated in detail. The influence of the acceleration and deceleration of the moving load on the dynamic response of the beam is also examined and discussed.

2. BFGSW Beam

A simply supported BFGSW beam with rectangular cross section (b × h) subjected to a load F₀, moving from left to right as depicted in Figure 1, is considered. The beam consists of three layers, a homogeneous core, and two FGM face layers with material properties varying in both the length and thickness directions. In the figure, the x-axis is chosen on the midplane, and z₀ = −h/2, z₁, z₂ and z₃ = h/2 are, respectively, the vertical coordinates of the bottom layer and the interfaces of the layers and the top layer.

[figure omitted; refer to PDF]

The beam is made from three materials, referred to as M1, M2, and M3. The core of the beam is pure M1, while the two face layers are bidirectional FGM with volume fraction of the constituents smoothly varying in both the x- and z-directions according to [43, 44]:$$\begin{array}{c}\text{(1)}& \left\{\begin{array}{c}{V}_1={\left(\frac{z-z_0}{z_1-z_0}\right)}^{n_z},\\ V_2=\left[1-{\left(\frac{z-z_0}{z_1-z_0}\right)}^{n_z}\right]\left[1-{\left(\frac{x}{L}\right)}^{n_x}\right],\hspace{1em}\mathrm{for}z\in \left[z_0,z_1\right],\\ V_3=\left[1-{\left(\frac{z-z_0}{z_1-z_0}\right)}^{n_z}\right]{\left(\frac{x}{L}\right)}^{n_x},\end{array}\right.V_1=1,V_2=V_3=0,\hspace{1em}\mathrm{for}z\in \left[z_1,z_2\right],\\ \left\{\begin{array}{c}{V}_1={\left(\frac{z-z_3}{z_2-z_3}\right)}^{n_z},\\ V_2=\left[1-{\left(\frac{z-z_3}{z_2-z_3}\right)}^{n_z}\right]\left[1-{\left(\frac{x}{L}\right)}^{n_x}\right],\hspace{1em}\mathrm{for}z\in \left[z_2,z_3\right],\\ V_3=\left[1-{\left(\frac{z-z_3}{z_2-z_3}\right)}^{n_z}\right]{\left(\frac{x}{L}\right)}^{n_x},\end{array}\right.\end{array}$$where L is the beam length; V₁, V₂, and V₃ are, respectively, the volume fraction of M1, M2, and M3; n_x and n_z are the axial and transverse grading indexes. It can be verified that (1) defines the volume fraction of the conventional unidirectional FGM Sandwich beams in [26, 27] if nx = 0 or M2 is identical to M3. Figure 2 shows the distribution in the thickness and longitudinal directions of V₁, V₂, and V₃ for n_x = n_z = 0.5 and z₂ = −z₁ = h/5.

[figures omitted; refer to PDF]

The effective material properties, P_f, such as the elastic moduli and mass density evaluated by Voigt’s model are of the following form [44]:$$\begin{array}{}\text{(2)}& P_f\left(x,z\right)=P_1{V}_1+P_2{V}_2+P_3{V}_3,\end{array}$$where P₁, P₂, and P₃ are the properties of the M1, M2, and M3, respectively.

Substituting (1) into (2), one gets$$\begin{array}{}\text{(3)}& P_f\left(x,z\right)=\left\{\begin{array}{cc}\left[P_1-P_{23}\left(x\right)\right]{\left(\frac{z-z_0}{z_1-z_0}\right)}^{n_z}+P_{23}\left(x\right),& \mathrm{for}z\in \left[z_0,z_1\right],\\ P_1,& \mathrm{for}z\in \left[z_1,z_2\right],\\ \left[P_1-P_{23}\left(x\right)\right]{\left(\frac{z-z_3}{z_2-z_3}\right)}^{n_z}+P_{23}\left(x\right),& \mathrm{for}z\in \left[z_2,z_3\right],\end{array}\right.\end{array}$$with$$\begin{array}{}\text{(4)}& P_{23}\left(x\right)=P_2-\left(P_2-P_3\right){\left(\frac{x}{L}\right)}^{n_x}.\end{array}$$

One can easily verify that if n_x = 0 or M2 is identical to M3, (3) returns to the effective properties of the unidirectional transverse FGM Sandwich beams made from M1 and M3. Figure 3 shows the variation of the effective Young’s modulus E_f and mass density ρ_f of the BFGSW beam made from a mixture of alumina (Al₂O₃) as M1, stainless steel (SUS304) as M2, and aluminum (Al) as M3 for n_x = n_z = 0.5 and n_x = n_z = 3. The properties of Al₂O₃, SUS304, and Al adopted from [32] are listed in Table 1.

[figures omitted; refer to PDF]

Table 1

Properties of constituent materials of BFGSW beam.

Materials	Role	E (GPa)	ρ (kg/m3)	υ
Alumina (Al2O3)	M1	380	3960	0.3
Steel (SUS304)	M2	210	7800	0.3
Aluminum (Al)	M3	70	2702	0.3

3. Mathematical Formulation

Based on the first-order shear deformation theory, the displacements in the x- and z-directions, u₁(x, z, t), and u₃(x, z, t), respectively, are given by$$\begin{array}{c}\text{(5)}& u_1\left(x,z,t\right)=u\left(x,t\right)-z\theta \left(x,t\right),u_3\left(x,z,t\right)=w\left(x,t\right),\end{array}$$where u(x, t) and $w$(x, t) are, respectively, the axial and transverse displacements of the point on the x-axis, θ(x, t) is the rotation of the cross section, and t is the time variable. Equation (5) leads to the axial train ε_xx and the shear strain γ_xz in the following forms:$$\begin{array}{c}\text{(6)}& {\epsilon }_{xx}=u_{,x}-z{\theta }_{,x},{\gamma }_{xz}=w_{,x}-\theta .\end{array}$$

In (6) and hereafter, a subscript comma is used to denote the derivative with respect to the variable which follows.

Constitutive equation based on linear behavior of the beam material is of the form$$\begin{array}{}\text{(7)}& \left\{\begin{array}{c}{\sigma }_{xx}\\ {\tau }_{xz}\end{array}\right\}=\left[\begin{array}{cc}{E}_f\left(x,z\right)& 0\\ 0& \psi G_f\left(x,z\right)\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\gamma }_{xz}\end{array}\right\},\end{array}$$where σ_xx and τ_xz are, respectively, the axial and shear stresses; E_f (x, z) and G_f (x, z) are the effective Young and shear moduli, defined by (3); $\psi$ is the shear correction factor, chosen by 5/6 for the beam with rectangular cross section herein.

From (6) and (7), one can write the strain energy of the beam in the form$$\begin{array}{c}\text{(8)}& U=\frac{1}{2}{\displaystyle {\int }_0^L{\displaystyle {\int }_A\left({\sigma }_{xx}{\epsilon }_{xx}+{\tau }_{xz}{\gamma }_{xz}\right)\text{d}A}\text{\hspace{0.17em}d}x}=\frac{1}{2}{\displaystyle {\int }_0^L\left[A_{11}{u}_{,x}^2-2A_{12}{u}_{,x}{\theta }_{,x}+A_{22}{\theta }_{,x}^2+\psi A_{33}{\left(w_{,x}-\theta \right)}^2\right]\text{d}x},\end{array}$$where A is the cross-sectional area; A₁₁, A₁₂, A₂₂, and A₃₃ are, respectively, the extensional, extensional-bending coupling, bending, and shear rigidities, defined as$$\begin{array}{c}\text{(9)}& \left(A_{11},A_{12},A_{22}\right)=b{\displaystyle {\int }_{-h/2}^{h/2}{E}_f\left(x,z\right)\left(1,z,z^2\right)\text{d}z}=b{\displaystyle \sum _{k=1}^3{\displaystyle {\int }_{z_{k-1}}^{z_k}{E}_f\left(x,z\right)\left(1,z,z^2\right)\text{d}z}},A_{33}=b{\displaystyle {\int }_{-h/2}^{h/2}{G}_f\left(x,z\right)\text{d}z}=b{\displaystyle \sum _{k=1}^3{\displaystyle {\int }_{z_{k-1}}^{z_k}{G}_f\left(x,z\right)\text{d}z}.}\end{array}$$

Substituting E_f (x, z) and G_f (x, z) from (3) into (9), one can write the rigidities in the form$$\begin{array}{}\text{(10)}& A_{ij}=A_{ij}^{\text{M}1}+A_{ij}^{\text{M}2}+A_{ij}^{\text{M}1\text{M}2}-A_{ij}^{\text{M}2\text{M}3}{\left(\frac{x}{L}\right)}^{n_x},\hspace{1em}\left(i,j=1,\dots ,3\right),\end{array}$$where $A_{ij}^{\text{M}1},A_{ij}^{\text{M}2},A_{ij}^{\text{M}1\text{M}2}$, and $A_{ij}^{\text{M}2\text{M}3}$ are, respectively, the rigidities of homogeneous and unidirectional transverse FGM beams made from M1, M2, and M3. These terms can be explicitly evaluated and their expressions are given by (A.1)–(A.4) in Appendix.

The kinetic energy resulted from (5) is of the form$$\begin{array}{c}\text{(11)}& T=\frac{1}{2}{\displaystyle {\int }_0^L{\displaystyle {\int }_A{\rho }_f\left(x,z\right)\left({\dot{u}}_1^2+{\dot{u}}_3^2\right)\text{d}A}\text{\hspace{0.17em}d}x}=\frac{1}{2}{\displaystyle {\int }_0^L\left[I_{11}\left({\dot{u}}^2+{\dot{w}}^2\right)-2I_{12}\dot{u}\dot{\theta }+I_{22}{\dot{\theta }}^2\right]\text{d}x},\end{array}$$where an over dot is used to denote the derivative with respect to the time variable t; ρ_f (x,z) is the effective mass density; I₁₁, I₁₂, and I₂₂ are the mass moments, defined as$$\begin{array}{}\text{(12)}& \left(I_{11},I_{12},I_{22}\right)=b{\displaystyle {\int }_{-h/2}^{h/2}{\rho }_f\left(x,z\right)\left(1,z,z^2\right)\text{d}z}=b{\displaystyle \sum _{k=1}^3{\displaystyle {\int }_{z_{k-1}}^{z_k}{\rho }_f\left(x,z\right)\left(1,z,z^2\right)\text{d}z}}.\end{array}$$

As the rigidities, the above mass moments can also be written in the form$$\begin{array}{}\text{(13)}& I_{ij}=I_{ij}^{\text{M}1}+I_{ij}^{\text{M}2}+I_{ij}^{\text{M}1\text{M}2}-I_{ij}^{\text{M}2\text{M}3}{\left(\frac{x}{L}\right)}^{n_x},\hspace{1em}\left(i,j=1,\dots ,3\right),\end{array}$$with $I_{ij}^{\text{M}1},I_{ij}^{\text{M}2},I_{ij}^{\text{M}1\text{M}2}$, and $I_{ij}^{\text{M}2\text{M}3}$ being given by (A.5)–(A.7) in Appendix.

The beam is assumed under action of the load F₀, moving from the left end to the right end as shown in Figure 1. The potential of the load F₀ is simply given by$$\begin{array}{}\text{(14)}& V=-{\displaystyle {\int }_0^L{F}_{0}w\left(x,t\right)\delta \left(x-s\left(t\right)\right)\text{d}x},\end{array}$$where $\delta \left(\cdot \right)$ is Dirac delta function; x is the abscissa, measured from the left end of the beam; s(t) is the function describing motion of the force F₀ at the time t, which can be calculated as$$\begin{array}{}\text{(15)}& s\left(t\right)=v_{0}t+\frac{a{t}^2}{2},\end{array}$$where $v_0$ is the initial speed of the moving load and a is the acceleration of the load which is assumed to be constant herein. The motion of the load F₀ is defined by the sign of the acceleration a; namely, if a = 0, the load is in a uniform motion (constant moving speed), a > 0 defines an acceleration motion, and a < 0 defines a deceleration motion.

4. Finite Element Formulation

The differential equation of motion for the beam can be obtained by applying Hamilton’s principle to (8), (11), and (14). However, due to the rigidities and mass moments, as seen from (10) and (13), functions of x, a closed-form solution for such equation is difficult to obtain. A finite element formulation is derived in this section for computing the dynamic response of the beam. To this end, the beam is assumed to be divided into a number of elements with length l. The vector of nodal displacements for a two-node beam element contains six components as$$\begin{array}{c}\text{(16)}& \mathbf{d}={\left\{\begin{array}{cc}{\mathbf{d}}_u& {\mathbf{d}}_b\end{array}\right\}}^T,\hspace{1em}\mathrm{with}{\mathbf{d}}_u={\left\{\begin{array}{cc}{u}_1& u_2\end{array}\right\}}^T,{\mathbf{d}}_b={\left\{\begin{array}{cccc}{w}_1& {\theta }_1& w_2& {\theta }_2\end{array}\right\}}^T,\end{array}$$where d_u and d_b denote the vectors of nodal axial and bending displacements, respectively; u_i, $w_i$, and θ_i (i = 1, 2) are, respectively, the axial displacement, transverse displacement, and rotation at the node i. The order of the nodal values is not necessary as in (16), but it is convenient to separate the nodal axial displacements from the bending ones. A superscript “T” in (16) and hereafter is used to denote the transpose of a vector or a matrix.

The displacements and rotation are interpolated from their nodal values according to$$\begin{array}{c}\text{(17)}& u={\mathbf{N}}_u{\mathbf{d}}_u,w={\mathbf{N}}_w{\mathbf{d}}_b,\\ \theta ={\mathbf{N}}_{\theta }{\mathbf{d}}_b,\end{array}$$where N_u = $\left\{\begin{array}{cc}{N}_{u1}& N_{u2}\end{array}\right\}$, ${\mathbf{N}}_w$ = $\left\{\begin{array}{cccc}{N}_{w1}& N_{w2}& N_{w3}& N_{w4}\end{array}\right\}$, and N_θ = $\left\{\begin{array}{cccc}{N}_{\theta 1}& N_{\theta 2}& N_{\theta 3}& N_{\theta 4}\end{array}\right\}$ are the matrices of interpolation functions for u, $w$, and θ, respectively. The following polynomials are employed as the interpolation functions in the present work.

For the axial displacement u,$$\begin{array}{c}\text{(18)}& N_{u1}=\frac{x}{l},N_{u2}=\frac{l-x}{l}.\end{array}$$

For the transverse displacement $w$,$$\begin{array}{c}\text{(19)}& N_{w1}=\frac{1}{\left(1+\varphi \right)}\left[2{\left(\frac{x}{l}\right)}^3-3{\left(\frac{x}{l}\right)}^2-\varphi \left(\frac{x}{l}\right)+\left(1+\varphi \right)\right],N_{w2}=\frac{l}{\left(1+\varphi \right)}\left[{\left(\frac{x}{l}\right)}^3-\left(2+\frac{\varphi }{2}\right){\left(\frac{x}{l}\right)}^2-\left(2+\frac{\varphi }{2}\right)\left(\frac{x}{l}\right)\right],\\ N_{w3}=-\frac{1}{\left(1+\varphi \right)}\left[2{\left(\frac{x}{l}\right)}^3-3{\left(\frac{x}{l}\right)}^2-\varphi \left(\frac{x}{l}\right)\right],\\ N_{w4}=\frac{l}{\left(1+\varphi \right)}\left[{\left(\frac{x}{l}\right)}^3-\left(1-\frac{\varphi }{2}\right){\left(\frac{x}{l}\right)}^2-\frac{\varphi }{2}\left(\frac{x}{l}\right)\right].\end{array}$$

For the rotation θ,$$\begin{array}{c}\text{(20)}& N_{\theta 1}=\frac{6}{\left(1+\varphi \right)l}\left[{\left(\frac{x}{l}\right)}^2-\left(\frac{x}{l}\right)\right],N_{\theta 2}=\frac{1}{\left(1+\varphi \right)}\left[3{\left(\frac{x}{l}\right)}^2-\left(4+\varphi \right)\left(\frac{x}{l}\right)+\left(1+\varphi \right)\right],\\ N_{\theta 3}=-\frac{6}{\left(1+\varphi \right)l}\left[{\left(\frac{x}{l}\right)}^2-\left(\frac{x}{l}\right)\right],\\ N_{\theta 4}=\frac{1}{\left(1+\varphi \right)}\left[3{\left(\frac{x}{l}\right)}^2-\left(2+\varphi \right)\left(\frac{x}{l}\right)\right].\end{array}$$

In (19) and (20), $\varphi =12A_{22}/\left(l^2\psi A_{33}\right)$ is the shear deformation parameter. The cubic and quadratic polynomials in (19) and (20) were firstly derived by Kosmatka in [45] for homogeneous beams, and they have been used by several authors [5, 39] to formulate the finite element formulations for analyzing FGM beams.

Using the interpolation functions (18)–(20), one can write the strain energy (8) in the form$$\begin{array}{}\text{(21)}& U=\frac{1}{2}{\displaystyle \sum ^{\text{NE}}{\mathbf{d}}^T\mathbf{k}\mathbf{d}},\end{array}$$where NE is the total number of elements used to discretize the beam and k is the element stiffness matrix, which can be written in submatrices as$$\begin{array}{}\text{(22)}& \mathbf{k}=\left[\begin{array}{cc}{\mathbf{k}}_{aa}& {\mathbf{k}}_{ab}\\ {\mathbf{k}}_{ab}^T& {\mathbf{k}}_{bb}+{\mathbf{k}}_{ss}\end{array}\right],\end{array}$$with k_aa, k_ab, k_bb, and k_ss being, respectively, the element stiffness matrices stemming from the axial stretching, axial-bending coupling, bending, and shear deformation. The expressions for these submatrices are as follows:$$\begin{array}{c}\text{(23)}& {\mathbf{k}}_{aa}={\displaystyle {\int }_0^l{\mathbf{N}}_{u,x}^T{A}_{11}{\mathbf{N}}_{u,x}\text{d}x},{\mathbf{k}}_{ab}=-{\displaystyle {\int }_0^l{\mathbf{N}}_{u,x}^T{A}_{12}{\mathbf{N}}_{\theta ,x}\text{d}x},\\ {\mathbf{k}}_{bb}={\displaystyle {\int }_0^l{\mathbf{N}}_{\theta ,x}^T{A}_{22}{\mathbf{N}}_{\theta ,x}\text{d}x},\\ {\mathbf{k}}_{ss}={\displaystyle {\int }_0^l{\left({\mathbf{N}}_{w,x}-{\mathbf{N}}_{\theta }\right)}^T\psi A_{33}\left({\mathbf{N}}_{w,x}-{\mathbf{N}}_{\theta }\right)\text{d}x}.\end{array}$$

Making use of (10), the integrals in (23) can be evaluated explicitly, and thus explicit expressions for the submatrices are obtained.

The kinetic energy (11) can also be written in the form$$\begin{array}{}\text{(24)}& T=\frac{1}{2}{\displaystyle \sum ^{\text{NE}}{\dot{\mathbf{d}}}^T\mathbf{m}\dot{\mathbf{d}}},\end{array}$$where the element mass matrix m can be written in submatrices as follows:$$\begin{array}{}\text{(25)}& \mathbf{m}=\left[\begin{array}{cc}{\mathbf{m}}_{uu}& {\mathbf{m}}_{u\theta }\\ {\mathbf{m}}_{u\theta }^T& {\mathbf{m}}_{ww}+{\mathbf{m}}_{\theta \theta }\end{array}\right],\end{array}$$with$$\begin{array}{c}\text{(26)}& {\mathbf{m}}_{uu}={\displaystyle {\int }_0^l{\mathbf{N}}_u^T{I}_{11}{\mathbf{N}}_u\text{d}x},{\mathbf{m}}_{ww}={\displaystyle {\int }_0^l{\mathbf{N}}_w^T{I}_{11}{\mathbf{N}}_w\text{d}x},\\ {\mathbf{m}}_{u\theta }=-{\displaystyle {\int }_0^l{\mathbf{N}}_u^T{I}_{12}{\mathbf{N}}_{\theta }\text{d}x},\\ {\mathbf{m}}_{\theta \theta }={\displaystyle {\int }_0^l{\mathbf{N}}_{\theta }^T{I}_{12}{\mathbf{N}}_{\theta }\text{d}x},\end{array}$$being, respectively, the element mass matrices resulting from the translations in axial and transverse directions, axial translation-sectional rotation coupling, and cross-sectional rotation. Making use of (13), explicit expressions for the submatrices in (26) can also be obtained.

The discrete equation of motion for the dynamic analysis of the beam can be written in the form [46]$$\begin{array}{}\text{(27)}& \mathbf{M}\ddot{\mathbf{D}}+\mathbf{K}\mathbf{D}={\mathbf{F}}_{\text{ex}},\end{array}$$where M and K are, respectively, the global mass and stiffness matrices, obtained by assembling the matrices m and k over the elements; D and $\ddot{\mathbf{D}}$ are, respectively, the vectors of nodal displacements and accelerations; F_ex is the vector of the nodal external force with the following form:$$\begin{array}{}\text{(28)}& {\mathbf{F}}_{\text{ex}}={\displaystyle \sum ^{\text{NE}}{\mathbf{f}}_{\text{ex}}},\hspace{1em}\mathrm{with}{\mathbf{f}}_{\text{ex}}={\left\{\begin{array}{cccccc}0& F_0{N}_{w1}& F_0{N}_{w2}& 0& F_0{N}_{w3}& F_0{N}_{w4}\end{array}\right\}}^T,\end{array}$$where f_ex is the element nodal force vector. Except for the element under the load F₀, the element nodal force vector f_ex is zero for all other elements, and the interpolation functions $N_{wi}\left(i=1,\dots ,4\right)$ in (28) are evaluated at the current position of the force F₀. Equation (27) can be solved by the Newmark method. The average acceleration method which ensures the numerically unconditional stability [46] is adopted herein.

5. Numerical Investigation

Numerical investigation is carried out in this section to examine the effects of various parameters such as the material grading indexes, layer thickness ratio, and moving load speed on the dynamic behavior of the beam. To this end, a simply supported beam with b = 0.5 m and h = 1 m, made from Al₂O₃, SUS304, and Al with the material data in Table 1, is employed herewith. Otherwise stated, an aspect ratio L/h = 20 is assumed. A uniform increment time step ∆t = ∆T/200 with ∆T being the total time, necessary for the load completely crossing the beam, is used for the Newmark procedure. For the convenience of discussion, the following dynamic magnification factor (DMF) is introduced:$$\begin{array}{}\text{(29)}& D_d=\max \left(\frac{w\left(\left(L/2\right),t\right)}{w_{\text{st}}}\right),\end{array}$$where $w_{\text{st}}$ = L³F₀/48E_sI is the static deflection of a fully SUS304 beam under the load F₀, acting at the midspan. Three numbers in brackets are used to denote the layer thickness ratio; for example, (2-1-2) means that the thickness ratio of the layers from bottom to top surfaces is 2 : 1 : 2.

5.1. Accuracy and Convergence Studies

Before computing the dynamic response of the beam, the accuracy and convergence of the derived element are firstly verified. Since there is no data on the present beam model in the literature, the verification is carried out for a unidirectional transverse FGM Sandwich beam made from Al₂O₃ and Al in [32, 33], a special case of the present BFGSW beam when M2 is identical to M3.

In Table 2, the frequency parameter, $\overline{\omega }=\left(\omega L^2\right)/\left(h\sqrt{{\rho }_{\text{Al}}/E_{\text{Al}}}\right)$ with ω being the fundamental frequency, obtained in the present work is compared to the result of Şimşek and Al-shujairi in [32], where a semianalytical method has been used. Very good agreement between the frequency parameter of the present work with that of [32] is noted from Table 2, regardless of the material index and the layer thickness ratio.

Table 2

Comparison of frequency parameter of unidirectional FGM Sandwich beam.

nz	Source	L/h = 5			L/h = 20
		(1-2-1)	(1-8-1)	(1-1-1)	(1-2-1)	(1-8-1)	(1-1-1)
0.5	Reference [32]	4.4711	4.8396	4.3192	4.6972	5.1064	4.5316
	Present work	4.4715	4.84	4.3196	4.6973	5.1065	4.5318
1	Reference [32]	4.1002	4.6762	3.8627	4.2881	4.9230	4.0318
	Present work	4.1006	4.6766	3.8632	4.2884	4.9231	4.0321
2	Reference [32]	3.7223	4.4510	3.4055	3.8760	4.7378	3.5379
	Present work	3.7228	4.5107	3.4061	3.8764	4.738	3.5383
5	Reference [32]	3.3652	4.3456	3.0039	3.4912	4.5550	3.1100
	Present work	3.3658	4.346	3.0045	3.4917	4.5552	3.1106
10	Reference [32]	3.2233	4.2727	2.8659	3.3397	4.4745	2.9651
	Present work	3.2239	4.2731	2.8666	3.3402	4.4747	2.9656

To verify the accuracy of the derived beam element in evaluating the dynamic response, Figure 4 compares the time histories for midspan dimensionless deflection of a transverse FGM Sandwich beam made from Al₂O₃ and Al obtained in the present paper with that of Songsuwan et al. [33] by using the Ritz method. Very good agreement between the present result with that of [33] is seen from Figure 4. Note that $w_{\text{st}}$ in Figure 4 is the maximum static deflection of the simply supported beam made from pure Al.

[figure omitted; refer to PDF]

Table 3 shows the convergence of the derived beam element in evaluating the fundamental frequency parameter of the BFGSW beam. As seen from the table, the convergence is achieved by using 26 elements, regardless of the material grading indexes and the layer thickness ratio. Note that the convergence of the results in Table 2 and also in Figure 4 has been achieved by using only 16 elements. Thus, the longitudinal variation of the material properties of the BFGSW beam makes the convergence of the element considerably slower. Because of this convergence result, 26 elements are used for all the computations reported in the following.

Table 3

Convergence of the element in evaluating frequencies of BFGSW beam (L/h = 20).

	nx	nz	NE = 16	NE = 18	NE = 20	NE = 22	NE = 24	NE = 26
(2-1-2)	0.5	0.3	4.7097	4.7096	4.7096	4.7096	4.7096	4.7096
		1	3.8637	3.8637	3.8636	3.8635	3.8634	3.8634
		3	3.2115	3.2113	3.2112	3.2111	3.2110	3.2110
	1	0.3	4.8366	4.8366	4.8366	4.8366	4.8365	4.8365
		1	3.8574	3.8573	3.8573	3.8572	3.8571	3.8571
		3	3.2521	3.2520	3.2519	3.2518	3.2517	3.2517
(2-2-1)	0.5	0.3	4.8576	4.8575	4.8575	4.8575	4.8575	4.8575
		1	4.1704	4.1703	4.1702	4.1702	4.1701	4.1701
		3	3.6098	3.6097	3.6096	3.6095	3.6094	3.6094
	1	0.3	3.6066	3.6066	3.6065	3.6095	3.6095	3.6095
		1	4.1599	4.15998	4.1598	4.1597	4.1597	4.1597
		3	3.6292	3.6290	3.6290	3.6289	3.6288	3.6288

5.2. Uniform Motion

The dynamic response of the BFGSW beam under uniform motion of the moving load (a = 0) is considered in this subsection. In Figure 5, the time histories for midspan dimensionless deflection of the BFGSW beam are depicted for various moving load speeds, different layer thickness ratios, and two pairs of the grading indexes, n_x = n_z = 0.5 and n_x = n_z = 3. In addition to the moving load speed, the figure also shows a significant influence of the grading indexes and the layer thickness ratio on the time histories of the beam. At a given moving load speed, the maximum deflection of the beam with grading indexes n_x = n_z = 3 is considerably higher than that of the beam with n_x = n_z = 0.5. The time at which the deflection attains the maximum value is also altered by the change of the grading indexes. The influence of the layer thickness ratio is also seen from Figure 5, where the maximum deflection is seen to be lower for the beam associated with a larger core thickness. This due to the fact that the beam with a larger core thickness has a higher percentage of Al₂O₃. Since Young’s modulus of Al₂O₃ is higher than that of the other constituents, the rigidities of the beam with a larger core thickness are higher, and this results in the lower deflection.

[figures omitted; refer to PDF]

In Figure 6, the transverse grading index n_z versus the DMF of the (2-2-1) beam is shown for different moving load speeds and two values of the axial index, n_x = 0.5 and n_x = 3. The influence of the transverse index n_z on the DMF of the BFGSW beam is similar to that of the transverse FGM sandwich beam, previously investigated in [33], and the DMF increases by the increase of the transverse index n_z, irrespective of the moving load speed. The increase of the DMF with the increase of the n_z can be explained by the decrease of the effective moduli, as can be seen from (3), and this leads to the decrease of the beam rigidities. The effect of the transverse index on the DMF, as seen from Figure 6, is more significant for n_z in range from 0 to 5. The axial index n_x and the moving load speed can change the amplitude of the DMF, but they hardly alter the variation of the factor with the transverse index n_z. Figure 6 also shows an important role of the moving load speed on the dynamic response, and the DMF of the beam increases with the increase of the moving load speed, regardless of the material grading indexes. The effects of the moving load speed and the material grading indexes on the DMF can be seen more clearly from Figures 7 and 8, where the variation of the DMF with the moving load speed and the material grading indexes are, respectively, depicted. As seen from the figures, the DMF increases with the increase of the transverse index n_z, but it decreases with the increase of the axial index n_x. The dependence of the DMF factor upon the indexes n_x and n_z of the BFGSW beam is similar to that of the bidirectional FGM beam as reported in [39]. The core thickness ratio, as seen from Figure 8, alters the amplitude of the DMF, but it hardly changes the dependence of this factor on the grading indexes. Based on the variation of the DMF with the grading indexes depicted in Figure 8, one can design the beam with a desired DMF by approximately choosing of the material grading indexes.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

In order to examine the influence of the aspect ratio on the dynamic response of the beam, Figure 9 shows the variation of the DMF with the moving load speed of the (2-1-2) and (2-2-1) beams for various values of the aspect ratio, L/h = 5, 10, 15, and 20, and for n_x = n_z = 2. The aspect ratio, as seen from the figure, has significant influence on the relation between the DMF and the moving load speed, and the DMF attains the maximum value at a higher moving load speed for the beam having a lower aspect ratio.

[figures omitted; refer to PDF]

5.3. Acceleration and Deceleration

The influence of the acceleration and deceleration of the moving load on the dynamic behavior of the BFGSW beam is examined herewith. It is assumed that for the deceleration motion the load enters the beam with a velocity $v$, and it exits the beam with a velocity $v$ = 0. On the other hand, for the acceleration motion, the load is assumed to enter the beam with a velocity $v$ = 0 and it exits the beam with a velocity $v$.

In Figure 10, the time histories for the midspan dimensionless deflection of the BFGSW beam under different motions of the moving load are shown for the (2-1-2) and (2-2-1) beams with n_x = n_z = 0.5 and $v$ = 100 m/s. The deceleration (a < 0) and acceleration (a > 0) motions, as seen from the figure, have a significant influence on the time history of the beam. As seen from Figure 10, the beam under the deceleration motion executives more vibration cycles than it does under the uniform and acceleration motions of the moving load. The maximum midspan deflection is also affected by the motion type, and the maximum deflection of the beam under the acceleration motion is considerably lower than that of the beam under the uniform and deceleration motions of the moving load. The effect of the deceleration and acceleration on the dynamic response of the BFGSW beam can also be seen from Figure 11, where the moving load speed versus the DMF is shown for the (2-1-2) beam with n_x = n_z = 0.5 and n_x = n_z = 3. For the most values of the moving load speed, the DMF of the beam under the acceleration motion of the moving load is considerably lower that than of the beam under the uniform and deceleration motions of the load.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

6. Conclusions

A BFGSW beam model made from three distinct materials has been proposed and its dynamic response to nonuniform motion of a moving load was studied. The beam consists of three layers, a homogeneous core, and two bidirectional FGM face sheets with material properties varying in both the thickness and length directions by the power gradation laws. The conventional transverse FGM Sandwich beam is a special case of the present BFGSW beam. Based on the first-order shear deformation beam theory, a finite beam element was formulated and employed in combination with Newmark method to compute the dynamic response of the beam. The element based on the quadratic and cubic interpolation for the rotation and transverse displacement is simple with the stiffness and mass matrices evaluated analytically. The accuracy of the formulated element has been confirmed through a comparison study. The obtained numerical results reveal that the material distribution and the layer thickness ratio play an important role on the dynamic behavior of the beam, and the beam can be designed to meet the desired DMF by choosing suitable grading indexes and a thickness ratio. It has been shown that the DMF increases with the increase of the transverse grading index, but it decreases with increasing the axial grading index. It was also shown that the dynamic characteristics, including the time histories for the midspan deflection and the DMF, are significantly influenced by the acceleration and deceleration of the moving load. The effects of the moving load speed, the material grading indexes, the layer thickness, and aspect ratios on the dynamic behavior of the beam have been examined in detail and highlighted. It is necessary to note that though the numerical investigation has been demonstrated for the simply supported beam only, the finite element formulation derived in the present work can be used to study dynamic response of BFGSW beam with other boundary conditions as well.

Acknowledgments

This work was supported by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under Grant no. 107.02-2018.23.

Appendix

Appendix

(1)

The rigidities $A_{ij}^{{\text{M}}_1},A_{ij}^{{\text{M}}_2},A_{ij}^{{\text{M}}_1{\text{M}}_2},A_{ij}^{{\text{M}}_2{\text{M}}_3}$ in (10) are as follows:

$$\begin{array}{c}\text{(A.1)}& A_{11}^{{\text{M}}_1}=\left(z_2-z_1\right)b{E}_1,A_{11}^{{\text{M}}_2}=\left(z_1-z_0+z_3-z_2\right)b{E}_2,\\ A_{11}^{{\text{M}}_1{\text{M}}_2}=\frac{\left(z_1-z_0+z_3-z_2\right)}{n_z+1}b{E}_{12},\\ A_{11}^{{\text{M}}_2{\text{M}}_3}=\frac{\left(z_1-z_0+z_3-z_2\right)n_z}{n_z+1}b{E}_{23},\\ \text{(A.2)}& A_{12}^{{\text{M}}_1}=\frac{z_2^2-z_1^2}{2}b{E}_1,A_{12}^{{\text{M}}_2}=\frac{z_1^2-z_0^2+z_3^2-z_2^2}{2}b{E}_2,\\ A_{12}^{{\text{M}}_1{\text{M}}_2}=\left[\frac{z_0\left(z_1-z_0\right)-z_3\left(z_2-z_3\right)}{n_z+1}+\frac{{\left(z_1-z_0\right)}^2-{\left(z_2-z_3\right)}^2}{n_z+2}\right]b{E}_{12},\\ A_{12}^{{\text{M}}_2{\text{M}}_3}=\left[\frac{z_1^2-z_0^2+z_3^2-z_2^2}{2}-\frac{z_0\left(z_1-z_0\right)-z_3\left(z_2-z_3\right)}{n_z+1}-\frac{{\left(z_1-z_0\right)}^2-{\left(z_2-z_3\right)}^2}{n_z+2}\right]b{E}_{23},\\ \text{(A.3)}& A_{22}^{{\text{M}}_1}=\frac{z_2^3-z_1^3}{3}b{E}_1,A_{22}^{{\text{M}}_2}=\frac{z_1^3-z_0^3+z_3^3-z_2^3}{3}b{E}_2,\\ A_{22}^{{\text{M}}_1{\text{M}}_2}=\left[\frac{z_0^2\left(z_1-z_0\right)-z_3^2\left(z_2-z_3\right)}{n_z+1}+2\frac{z_0{\left(z_1-z_0\right)}^2-z_3{\left(z_2-z_3\right)}^2}{n_z+2}\right.+\left.\frac{{\left(z_1-z_0\right)}^3-{\left(z_2-z_3\right)}^3}{n_z+3}\right]b{E}_{12},\\ A_{22}^{{\text{M}}_2{\text{M}}_3}=\left[\frac{z_1^3-z_0^3+z_3^3-z_2^3}{3}-\frac{z_0^2\left(z_1-z_0\right)-z_3^2\left(z_2-z_3\right)}{n_z+1}\right.-2\left.\frac{z_0{\left(z_1-z_0\right)}^2-z_3{\left(z_2-z_3\right)}^2}{n_z+2}-\frac{{\left(z_1-z_0\right)}^3-{\left(z_2-z_3\right)}^3}{n_z+3}\right]b{E}_{23},\\ \text{(A.4)}& A_{33}^{{\text{M}}_1}=\left(z_2-z_1\right)b{G}_1,A_{33}^{{\text{M}}_2}=\left(z_1-z_0+z_3-z_2\right)b{G}_2,\\ A_{33}^{{\text{M}}_1{\text{M}}_2}=\frac{\left(z_1-z_0+z_3-z_2\right)}{n_z+1}b{G}_{12},\\ A_{33}^{{\text{M}}_2{\text{M}}_3}=\frac{\left(z_1-z_0+z_3-z_2\right)}{n_z+1}b{G}_{23},\end{array}$$with $E_{12}=E_1-E_2,E_{23}=E_2-E_3,G_{12}=G_1-G_2,G_{23}=G_2-G_3$.

(2)

Mass moments $I_{ij}^{{\text{M}}_1},I_{ij}^{{\text{M}}_2},I_{ij}^{{\text{M}}_1{\text{M}}_2},I_{ij}^{{\text{M}}_2{\text{M}}_3}$ in (13) are as follows:

$$\begin{array}{c}\text{(A.5)}& I_{11}^{{\text{M}}_1}=\left(z_2-z_1\right)b{\rho }_1,I_{11}^{{\text{M}}_2}=\left(z_1-z_0+z_3-z_2\right)b{\rho }_2,\\ I_{11}^{{\text{M}}_1{\text{M}}_2}=\frac{\left(z_1-z_0+z_3-z_2\right)}{n_z+1}b{\rho }_{12},\\ I_{11}^{{\text{M}}_2{\text{M}}_3}=\frac{\left(z_1-z_0+z_3-z_2\right)n_z}{n_z+1}b{\rho }_{23},\\ \text{(A.6)}& I_{12}^{{\text{M}}_1}=\frac{z_2^2-z_1^2}{2}b{\rho }_1,I_{12}^{{\text{M}}_2}=\frac{z_1^2-z_0^2+z_3^2-z_2^2}{2}b{\rho }_2,\\ I_{12}^{{\text{M}}_1{\text{M}}_2}=\left[\frac{z_0\left(z_1-z_0\right)-z_3\left(z_2-z_3\right)}{n_z+1}+\frac{{\left(z_1-z_0\right)}^2-{\left(z_2-z_3\right)}^2}{n_z+2}\right]b{\rho }_{12},\\ I_{12}^{{\text{M}}_2{\text{M}}_3}=\left[\frac{z_1^2-z_0^2+z_3^2-z_2^2}{2}-\frac{z_0\left(z_1-z_0\right)-z_3\left(z_2-z_3\right)}{n_z+1}-\frac{{\left(z_1-z_0\right)}^2-{\left(z_2-z_3\right)}^2}{n_z+2}\right]b{\rho }_{23},\\ \text{(A.7)}& I_{22}^{{\text{M}}_1}=\frac{z_2^3-z_1^3}{3}b{\rho }_1,I_{22}^{{\text{M}}_2}=\frac{z_1^3-z_0^3+z_3^3-z_2^3}{3}b{\rho }_2,\\ I_{22}^{{\text{M}}_1{\text{M}}_2}=\left[\frac{z_0^2\left(z_1-z_0\right)-z_3^2\left(z_2-z_3\right)}{n_z+1}+2\frac{z_0{\left(z_1-z_0\right)}^2-z_3{\left(z_2-z_3\right)}^2}{n_z+2}\right.+\left.\frac{{\left(z_1-z_0\right)}^3-{\left(z_2-z_3\right)}^3}{n_z+3}\right]b{\rho }_{12},\\ I_{22}^{{\text{M}}_2{\text{M}}_3}=\left[\frac{z_1^3-z_0^3+z_3^3-z_2^3}{3}-\frac{z_0^2\left(z_1-z_0\right)-z_3^2\left(z_2-z_3\right)}{n_z+1}\right.-2\left.\frac{z_0{\left(z_1-z_0\right)}^2-z_3{\left(z_2-z_3\right)}^2}{n_z+2}-\frac{{\left(z_1-z_0\right)}^3-{\left(z_2-z_3\right)}^3}{n_z+3}\right]b{\rho }_{23},\end{array}$$with ${\rho }_{12}={\rho }_1-{\rho }_2,{\rho }_{23}={\rho }_2-{\rho }_3$.

References

[1] M. Koizumi, "FGM activities in Japan," Composites Part B: Engineering, vol. 28 no. 1-2,DOI: 10.1016/s1359-8368(96)00016-9, 1997.

[2] S. A. Sina, H. M. Navazi, H. Haddadpour, "An analytical method for free vibration analysis of functionally graded beams," Materials and Design, vol. 30 no. 3, pp. 741-747, DOI: 10.1016/j.matdes.2008.05.015, 2009.

[3] A. Mahi, E. A. Adda Bedia, A. Tounsi, I. Mechab, "An analytical method for temperature-dependent free vibration analysis of functionally graded beams with general boundary conditions," Composite Structures, vol. 92 no. 8, pp. 1877-1887, DOI: 10.1016/j.compstruct.2010.01.010, 2010.

[4] A. E. Alshorbagy, M. A. Eltaher, F. F. Mahmoud, "Free vibration characteristics of a functionally graded beam by finite element method," Applied Mathematical Modelling, vol. 35 no. 1, pp. 412-425, DOI: 10.1016/j.apm.2010.07.006, 2011.

[5] A. Shahba, R. Attarnejad, M. T. Marvi, 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, vol. 42 no. 4, pp. 801-808, DOI: 10.1016/j.compositesb.2011.01.017, 2011.

[6] N. Wattanasakulpong, B. Gangadhara Prusty, D. W. Kelly, M. Hoffman, "Free vibration analysis of layered functionally graded beams with experimental validation," Materials & Design (1980–2015), vol. 36, pp. 182-190, DOI: 10.1016/j.matdes.2011.10.049, 2012.

[7] M. Şimşek, T. Kocatürk, "Free and forced vibration of a functionally graded beam subjected to a concentrated moving harmonic load," Composite Structures, vol. 90 no. 4, pp. 465-473, DOI: 10.1016/j.compstruct.2009.04.024, 2009.

[8] M. Şimşek, "Vibration analysis of a functionally graded beam under a moving mass by using different beam theories," Composite Structures, vol. 92 no. 4, pp. 904-917, DOI: 10.1016/j.compstruct.2009.09.030, 2010.

[9] M. Şimşek, T. Kocatürk, Ş. D. Akbas, "Dynamic behavior of an axially functionally graded beam under action of a moving harmonic load," Composite Structures, vol. 94 no. 8, pp. 2358-2364, DOI: 10.1016/j.compstruct.2012.03.020, 2012.

[10] S. M. R. Khalili, A. A. Jafari, S. A. Eftekhari, "A mixed Ritz-DQ method for forced vibration of functionally graded beams carrying moving loads," Composite Structures, vol. 92 no. 10, pp. 2497-2511, DOI: 10.1016/j.compstruct.2010.02.012, 2010.

[11] K. Rajabi, M. H. Kargarnovin, M. Gharini, "Dynamic analysis of a functionally graded simply supported Euler-Bernoulli beam subjected to a moving oscillator," Acta Mechanica, vol. 224 no. 2, pp. 425-446, DOI: 10.1007/s00707-012-0769-y, 2013.

[12] B. S. Gan, T.-H. Trinh, T.-H. Le, D.-K. Nguyen, "Dynamic response of non-uniform Timoshenko beams made of axially FGM subjected to multiple moving point loads," Structural Engineering and Mechanics, vol. 53 no. 5, pp. 981-995, DOI: 10.12989/sem.2015.53.5.981, 2015.

[13] Y. Wang, D. Wu, "Thermal effect on the dynamic response of axially functionally graded beam subjected to a moving harmonic load," Acta Astronautica, vol. 127, pp. 171-181, DOI: 10.1016/j.actaastro.2016.05.030, 2016.

[14] D. K. Nguyen, V. T. Bui, "Dynamic analysis of functionally graded Timoshenko beams in thermal environment using a higher-order hierarchical beam element," Mathematical Problems in Engineering, vol. 2017,DOI: 10.1155/2017/7025750, 2017.

[15] I. Esen, M. A. Koc, Y. C ay, "Finite element formulation and analysis of a functionally graded Timoshenko beam subjected to an accelerating mass including inertial effects of the mass," Latin American Journal of Solids and Structures, vol. 15 no. 10,DOI: 10.1590/1679-78255102, 2018.

[16] Y. Fukui, "Fundamental investigation of functionally gradient material manufacturing system using centrifugal force," JSME International Journal. Ser. 3, Vibration, Control Engineering, Engineering for Industry, vol. 34 no. 1, pp. 144-148, DOI: 10.1299/jsmec1988.34.144, 1991.

[17] A. Chakraborty, S. Gopalakrishnan, J. N. Reddy, "A new beam finite element for the analysis of functionally graded materials," International Journal of Mechanical Sciences, vol. 45 no. 3, pp. 519-539, DOI: 10.1016/s0020-7403(03)00058-4, 2003.

[18] N. A. Apetre, B. V. Sankar, D. R. Ambur, "Analytical modeling of sandwich beams with functionally graded core," Journal of Sandwich Structures & Materials, vol. 10 no. 1, pp. 53-74, DOI: 10.1177/1099636207081111, 2008.

[19] O. Rahmani, S. M. R. Khalili, K. Malekzadeh, H. Hadavinia, "Free vibration analysis of sandwich structures with a flexible functionally graded syntactic core," Composite Structures, vol. 91 no. 2, pp. 229-235, DOI: 10.1016/j.compstruct.2009.05.007, 2009.

[20] S. C. Pradhan, T. Murmu, "Thermo-mechanical vibration of an FGM sandwich beam under variable elastic foundations using differential quadrature method," Journal of Sound and Vibration, vol. 321 no. 1-2, pp. 342-362, DOI: 10.1016/j.jsv.2008.09.018, 2009.

[21] A. M. Zenkour, M. N. M. Allam, M. Sobhy, "Bending analysis of FG viscoelastic sandwich beams with elastic cores resting on Pasternak’s elastic foundations," Acta Mechanica, vol. 212 no. 3-4, pp. 233-252, DOI: 10.1007/s00707-009-0252-6, 2010.

[22] Z. Su, G. Jin, Y. Wang, X. Ye, "A general Fourier formulation for vibration analysis of functionally graded sandwich beams with arbitrary boundary condition and resting on elastic foundations," Acta Mechanica, vol. 227 no. 5, pp. 1493-1514, DOI: 10.1007/s00707-016-1575-8, 2016.

[23] M. C. Amirani, S. M. R Khalili, N. Nemati, "Free vibration analysis of sandwich beam with FG core using the element free Galerkin method," Composite Structures, vol. 90 no. 3, pp. 373-379, DOI: 10.1016/j.compstruct.2009.03.023, 2009.

[24] Y. Yang, C. C. Lam, K. P. Kou, V. P. Iu, "Free vibration analysis of the functionally graded sandwich beams by a meshfree boundary-domain integral equation method," Composite Structures, vol. 117, pp. 32-39, DOI: 10.1016/j.compstruct.2014.06.016, 2014.

[25] T. P. Vo, H.-T. Thai, T.-K. Nguyen, A. Maheri, J. Lee, "Finite element model for vibration and buckling of functionally graded sandwich beams based on a refined shear deformation theory," Engineering Structures, vol. 64, pp. 12-22, DOI: 10.1016/j.engstruct.2014.01.029, 2014.

[26] T. P. Vo, H.-T. Thai, T.-K. Nguyen, F. Inam, J. Lee, "A quasi-3D theory for vibration and buckling of functionally graded sandwich beams," Composite Structures, vol. 119,DOI: 10.1016/j.compstruct.2014.08.006, 2015.

[27] T. P. Vo, H.-T. Thai, T.-K. Nguyen, F. Inam, J. Lee, "Static behaviour of functionally graded sandwich beams using a quasi-3D theory," Composites Part B: Engineering, vol. 68, pp. 59-74, DOI: 10.1016/j.compositesb.2014.08.030, 2015.

[28] T.-K. Nguyen, T. Truong-Phong Nguyen, T. P. Vo, H.-T. Thai, "Vibration and buckling analysis of functionally graded sandwich beams by a new higher-order shear deformation theory," Composites Part B: Engineering, vol. 76, pp. 273-285, DOI: 10.1016/j.compositesb.2015.02.032, 2015.

[29] T.-K. Nguyen, B.-D. Nguyen, "A new higher-order shear deformation theory for static, buckling and free vibration analysis of functionally graded sandwich beams," Journal of Sandwich Structures & Materials, vol. 17 no. 6, pp. 613-631, DOI: 10.1177/1099636215589237, 2015.

[30] T.-K. Nguyen, T. P. Vo, B.-D. Nguyen, J. Lee, "An analytical solution for buckling and vibration analysis of functionally graded sandwich beams using a quasi-3D shear deformation theory," Composite Structures, vol. 156, pp. 238-252, DOI: 10.1016/j.compstruct.2015.11.074, 2016.

[31] R. Bennai, H. A. Atmane, A. Tounsi, "A new higher-order shear and normal deformation theory for functionally graded sandwich beams," Steel and Composite Structures, vol. 19 no. 3, pp. 521-546, DOI: 10.12989/scs.2015.19.3.521, 2015.

[32] M. Şimşek, M. Al-shujairi, "Static, free and forced vibration of functionally graded (FG) sandwich beams excited by two successive moving harmonic loads," Composites Part B, vol. 108, pp. 18-34, 2017.

[33] W. Songsuwan, M. Pimsarn, N. Wattanasakulpong, "Dynamic responses of functionally graded sandwich beams resting on elastic foundation under harmonic moving loads," International Journal of Structural Stability and Dynamics, vol. 18 no. 9,DOI: 10.1142/S0219455418501122, 2018.

[34] M. Şimşek, "Bi-directional functionally graded materials (BFGMs) for free and forced vibration of Timoshenko beams with various boundary conditions," Composite Structures, vol. 133, pp. 969-997, DOI: 10.1016/j.compstruct.2015.08.021, 2015.

[35] D. Hao, C. Wei, "Dynamic characteristics analysis ofbi-directional functionally graded Timoshenko beams," Composite Structures, vol. 141, pp. 253-263, DOI: 10.1016/j.compstruct.2016.01.051, 2016.

[36] Z.-h. Wang, X.-h. Wang, G.-d. Xu, S. Cheng, T. Zeng, "Free vibration of two-directional functionally graded beams," Composite Structures, vol. 135, pp. 191-198, DOI: 10.1016/j.compstruct.2015.09.013, 2016.

[37] T. A. Huynh, X. Q. Lieu, J. Lee, "NURBS-based modeling of bidirectional functionally graded Timoshenko beams for free vibration problem," Composite Structures, vol. 160, pp. 1178-1190, DOI: 10.1016/j.compstruct.2016.10.076, 2017.

[38] A. Karamanli, "Bending behaviour of two directional functionally graded sandwich beams by using a quasi-3d shear deformation theory," Composite Structures, vol. 174, pp. 70-86, 2017.

[39] D. K. Nguyen, Q. H. Nguyen, T. T. Tran, V. T. Bui, "Vibration of bi-dimensional functionally graded Timoshenko beams excited by a moving load," Acta Mechanica, vol. 228 no. 1, pp. 141-155, DOI: 10.1007/s00707-016-1705-3, 2017.

[40] D. K. Nguyen, T. T. Tran, "Free vibration of tapered BFGM beams using an efficient shear deformable finite element model," Steel and Composite Structures, vol. 29 no. 3, pp. 363-377, 2018.

[41] A. A. Bousahla, F. Bourada, S. R. Mahmoud, "Buckling and dynamic behavior of the simply supported CNT-RC beams using an integral-first shear deformation theory," Computers and Concrete, vol. 25 no. 2, pp. 155-166, 2020.

[42] H. Matouk, A. A. Bousahla, H. Heireche, "Investigation on hygro-thermal vibration of P-FG and symmetric S-FG nanobeam using integral Timoshenko beam theory," Advances in Nano Research, vol. 8 no. 4, pp. 293-305, 2020.

[43] M. Nemat-Alla, K. I. E. Ahmed, I. Hassab-Allah, "Elastic-plastic analysis of twodimensional functionally graded materials under thermal loading," International Journal of Solids and Structures, vol. 46 no. 14-15, pp. 2774-2786, DOI: 10.1016/j.ijsolstr.2009.03.008, 2009.

[44] T. V. Do, D. K. Nguyen, N. D. Duc, D. H. Doan, T. Q. Bui, "Analysis of bi-directional functionally graded plates by FEM and a new third-order shear deformation plate theory," Thin-Walled Structures, vol. 119, pp. 687-699, 2017.

[45] J. B. Kosmatka, "An improved two-node finite element for stability and natural frequencies of axial-loaded Timoshenko beams," Computers & Structures, vol. 57 no. 1, pp. 141-149, DOI: 10.1016/0045-7949(94)00595-t, 1995.

[46] R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, Concepts and Applications of Finite Element Analysis, 2002.