Dynamic characteristics analysis of bi-directional functionally graded Timoshenko beams

Accepted Manuscript Dynamic characteristics analysis of bi-directional functionally graded Timoshenko beams Deng Hao, Cheng Wei PII: DOI: Reference:

S0263-8223(16)00064-7 http://dx.doi.org/10.1016/j.compstruct.2016.01.051 COST 7156

To appear in:

Composite Structures

Please cite this article as: Hao, D., Wei, C., Dynamic characteristics analysis of bi-directional functionally graded Timoshenko beams, Composite Structures (2016), doi: http://dx.doi.org/10.1016/j.compstruct.2016.01.051

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Dynamic characteristics analysis of bi-directional functionally graded Timoshenko beams DengHao a * , ChengWei a a School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China;

Abstract: In this paper, the motion differential equations of the bi-directional functionally graded Timoshenko beam are established using Hamilton’s principle. The material properties of the beam change exponentially in both axial and thickness directions. Using the variable substitution method, state space differential equations of the structure are established. First, the dynamic stiffness matrix is formed using the traditional method. Then, the author proposes a new method to directly form an exact dynamic stiffness matrix by using state space differential equations, and this method is compared with the traditional dynamic stiffness matrix method. At the same time, the natural frequency of the structure is computed by combining the Wittrick-William algorithm with a non-iterative algorithm. The influence of gradient parameters α , β on the fundamental frequency, mode shape and frequency response function is analysed through the establishment of the dynamic stiffness matrix of the overall structure. Finally, using the Lagrange equation and the method of modal superposition, structural dynamic differential equations under a harmonic moving load are derived. Using the precise integration method, the dynamic response of the displacement is computed and the influence of gradient parameters α , β on the dynamic response is analysed. Key words: vibration; bi-directional functionally graded material; dynamic characteristics; exact solution; Timoshenko beam

1 Introduction The properties of functionally graded materials can be varied continuously in space, so these materials can effectively eliminate the problem of material parameters mismatching at interfaces. Because a functionally graded material can improve the strength and stiffness of a structure, and also has a good resistance to high temperature, these advantages lead it to be more widely used than traditional isotropic homogeneous materials. The gradient property of this material is suitable for industrial design with special requirements, so functionally graded materials are widely used in aerospace, nuclear energy, biomedical engineering and other fields [1-2]. There have been many advanced mechanics studies concerning the design, manufacture and service process of functionally graded materials. Inhomogeneity of material properties makes mechanical analysis very difficult. Existing mechanical concepts, theories, calculation methods and test methods for homogeneous materials are no longer suitable for functionally graded materials. Therefore, these areas need to be explored and developed. The further study of these mechanical problems is a prerequisite for the further application of functionally graded materials. At the same time, it also promotes the development of research on the mechanics of the inhomogeneous medium. Thus, it is of great theoretical significance to study the key mechanical problems of functionally graded materials and structures. Material properties varying towards the thickness direction have been intensively studied in recent years [3-6]. Sankar [7] studied the static deformation of a functionally graded beam subjected to transverse loads under simply supported boundary conditions. The dynamic response of a functionally graded beam with various end conditions was studied by Li [8]. *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

Simsek and Kocaturk [9] studied the dynamic characteristics of functionally graded beams with simply supported ends under concentrated moving harmonic loads. Chakraverty [10] used the Rayleigh-Ritz method to obtain the free vibration frequencies of functionally graded Timoshenko beams, and the results were compared with Euler- Bernoulli beams. For material properties graded along the axial direction, Elishakoff [11] proposed a semi-inverse method to obtain the exact solution based on the Euler-Bernoulli beam theory. Analytical solutions of Euler-Bernoulli beams with exponential bending stiffness and mass density under different boundary conditions have been obtained in [12-13]. Wang [14-15] computed the resonance frequencies of an Euler-Bernoulli beam with constant thickness and exponentially decaying width carrying a tip mass. Axial functionally graded tapered beams were analysed by Rajasekaran [16] using differential transformation and differential quadrature methods. Esmailzadeh [17] applied the Frobenius series method to analyse the vibration and stability of a Timoshenko beam subjected to an axial load. Huang [18] et al. used a new method to compute natural frequencies of axially functionally graded Timoshenko beams with non-uniform cross-section. Tang [19] obtained analytical solutions for a functionally graded Timoshenko beam with exponentially decaying material properties under different boundary conditions. In recent years, modified couple stress theory and various shear deformation beam theories have been proposed. M.Simsek [20] applied different high-order beam theories to analyse the fundamental frequency of functionally graded beams. Salamat-talab M [21] analysed the static and dynamic characteristics of a third-order shear deformation FG micro beam based on modified couple stress theory. Vo TP [22] established a finite element model for vibration and buckling of functionally graded sandwich beams based on a refined shear deformation theory. Zhang B [23] established the size-dependent functionally graded beam model based on an improved third-order shear deformation theory. Simsek M [24] investigated the buckling of a functionally graded microbeam embedded in elastic medium by using modified couple stress theories and a unified higher order beam theory. Ke LL [25] analysed the size effect on the dynamic stability of functionally graded microbeams based on a modified couple stress theory. In advanced machines, such as modern aerospace shuttles and craft, the temperature field develops in two or three directions, which means that conventional functionally graded materials may not be useful in the design of such structures [26]. As a result of this demand, the material properties of functionally graded materials are required to be graded in two or three directions. However, a large number of studies on functionally graded materials have been confined to the variation of material properties along the thickness direction or the axial direction, but studies on bi-directional functionally graded Timoshenko beams are very rare. Lü [27] proposed a novel elasticity solution combining the state-space method with differential quadrature method to solve static deformations of bi-directional functionally graded beams whose properties are varied exponentially. This semi-analytical and semi-numerical method was improved by Nie et al. [28-29], who proposed that the first derivative of displacement and displacement can be used as the state variable. This method deals more easily with certain boundary conditions than the state space method, which mixes stress and displacement as the state variable. This method has been applied to solve bi-directional graded beams [30] and two-dimensional bi-directional functionally graded circular plates [31]. Recently, Mesut Simsek [32] applied the Rayleigh method and Lagrange equations to establish the motion differential equations of a bi-directional functionally graded Timoshenko (BFGT) beam under the action of moving loads, and the dynamic response of the displacement under various boundary conditions *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

was computed using the polynomial trial functions that satisfied the boundary conditions. However, this method is a numerical method and its accuracy is not high. In this paper, to obtain the exact solution of a bi-directional functionally graded Timoshenko beam, the authors propose to combine the state space method with the dynamic stiffness method to solve motion differential equations. The concept of the dynamic stiffness matrix method [33-39] was first proposed by Kolousek [40] in the early 1940s. This method is a powerful tool to solve the problem of structural vibration in engineering, especially when higher order natural frequencies and higher accuracy are needed. The dynamic stiffness matrix method is also commonly referred to as an exact method. This is because, unlike the traditional finite element method or other approximate methods, it can obtain infinite natural frequencies and mode shapes of arbitrary precision. This method only needs to establish the initial displacement field assumption, and can be used to solve the motion differential equation accurately without any need to introduce other hypotheses or approximations, so it can give accurate results without regard to the number of elements. The dynamic stiffness matrix has the properties of element mass and element stiffness, and the matrix element is the transcendental function of frequency. The matrix is made up of analytical solutions of the motion differential equations, and this problem is often referred to as the transcendental eigenvalue problem. H.Su [41] formed the dynamic stiffness matrix of the functionally graded Timoshenko beam element using the traditional dynamic stiffness method based on their earlier research, and the influence of gradient parameters on the natural frequency of the structure was also analysed. Previously, however, the exact dynamic stiffness matrix of the bi-functionally graded Timoshenko beam whose properties are varied exponentially has not been obtained. High-order beam theories are able to represent the section warping in the deformed configuration and have higher accuracy for high order modes, while acquiring the exact solutions is complex and difficult. For simplicity, Timoshenko beam theory is considered in this paper. Timoshenko beam theory requires shear correction factors and is suitable for composite beams for which the shear correction factors can be determined. On the basis of previous studies, the authors establish the exact element dynamic stiffness matrix of the bi-directional functionally graded Timoshenko beam. In addition, the influence of the gradient parameters α , β on the fundamental frequency, mode shape and dynamic response under the action of a harmonic moving load is analysed.

2. Theory and formulation 2.1 Motion differential equations A bi-directional functionally graded Timoshenko beam (BFGT) with a rectangular cross section is shown in Fig. 1; the beam has a length L , width b and thickness h . The material properties of the beam are elastic modulus E , shear modulus G , Poisson’s ratio ν and mass density ρ . The material properties at positions 1,2,3,4 are Plb , Plt , Prb , Prt . The material properties obey an exponential distribution along the axial direction and the thickness direction, except for Poisson's ratio, as shown in Fig. 1. The material gradient parameter along the axis is β ,and the gradient parameter along the thickness is α : z 1 α ⋅ + 

Plt = e

h 2

⋅ Plb , Prb = e

β⋅

y L

⋅ Plb

(1)

Therefore, the material properties of the BFGT can be expressed as z 1 y α  +  + β

P( y , z ) = Plb e

h 2

L

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

(2)

z 4, Prt

2, Plt

y

3, Prb

1, Plb

Fig.1 Bi-directional functionally graded Timoshenko beam

Consider the deformation in the y-z plane of the functionally graded Timoshenko beam. The displacement of points on the neutral axis in the y direction is v and in the z direction is w . The displacements of any point on the cross section perpendicular to the neutral axis in the y direction and z direction are u y and u z , respectively. Thus, u y ( y , z , t ) = v ( y , t ) − z ⋅ w, y ( y , t ) + f ( z ) ⋅ψ ( y, t ) uz ( y, z, t ) = w( y, t )

ψ ( y , t ) = w, y ( y , t ) − φ ( y , t )

（3）

where φ ( y , t ) is the cross section rotation angle and f ( z ) is the distribution function of shear stress and shear strain of the cross section in the thickness direction [20]. For an Euler-Bernoulli beam: f ( z ) = 0; （4） For a Timoshenko beam： f ( z ) = z;

（5）

According to Timoshenko beam theory, substituting Eq.(5) into Eq.(3) yields: u y ( y , z , t ) = v ( y , t ) − z ⋅ φ ( y, t )

（6）

Strain and displacement should satisfy the following relation according to elasticity theory.  ∂u y   ∂v ∂φ    ∂y − z ⋅ ∂y  ε yy   ∂y     = =  γ ∂ u ∂ w ∂ u yz    y + z  −φ   ∂y   ∂y  ∂z

（7）

where ε yy is normal strain in the y direction and γ yz is the shear strain. According to Hooke’s law, the stress and strain of any point at the beam cross section should satisfy: σ yy   E ( y , z ) 0  ε yy  τ  =    G ( y , z )  γ yz   yz   0

（8）

Using the above relationships, the potential energy U and kinetic energy T of the bi-directional functionally graded beam are: 1 L (σ yyε yy +τ yzγ yz )dAdy 2 ∫0 ∫A y y y y y y  β β β β β 1 L β L v φ + A e Lφ2 + A e L w2 − 2 A e L w φ + A e Lφ 2 dy = ∫  A0e L v,2y − 2 Ae  1 , y , y 2 , y 3 , y 3 , y 3 2 0  

U=

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

（9）

 ∂u y 2  ∂u  2  1 L z  ρ ( y , z )  +  dAdy 2 ∫0 ∫A  ∂t   ∂t   y y y β β 1 L  β  ∂w ∂v  ∂v ∂φ ∂φ  = ∫  I 0 e L  ( ) 2 + ( )2  − 2 I1e L ( ) ⋅ ( ) + I 2 e L ⋅ ( )2 dy 2 0  ∂t  ∂t ∂t ∂t   ∂t

T=

（10）

(), y represents the first order derivative with respect to the coordinate y . The parameters

Ii (i = 0,1, 2) and A j ( j = 0,1, 2, 3) are defined as: z 1 α ⋅ + 

Ii = ∫ zi e

 h 2 ρ

z 1 α ⋅ + 

i lb dA, Ai = ∫ z e

 h 2  E dA( i lb

= 0,1, 2)

(11)

z 1 α ⋅ + 

A3 = ∫ κ e

 h 2  G dA lb

where κ is the shear coefficient whose value here is κ = 5 / 6 . Substituting for the variables w, v, φ : W0 = w ⋅ e

β

y L ,V 0

= v ⋅e

β

y L ,Φ

= φ ⋅e

0

β

y L

(12)

Thus, w′e

φ ′e

β

y L

= W0′ −

y L

β

= Φ ′0 −

β L

β L

W0 , v ′e

β

y L

= V0′ −

β L

V0

（13）

Φ0

Hamilton’s principle states: t2

δ ∫ (T − U )dt = 0

(14)

t1

where t1 and t2 are the time interval and δ is the usual variational operator, respectively. The motion differential equations of the BFGT can be obtained by substituting Eqs.(9)-(11) into Eq.(14) and noting that δ u, δ v, δφ are completely arbitrary, β

−I0

eL

β

I1e L

y

β

β

β y

β

y

y

∂ ( e v ′) ∂ ( e φ ′) v

+ I1φ

e L + A0 − A1 =0 ∂y ∂y y

L

y

β

L

β

β

y

β

y y y ∂ ( e L v ′) ∂ ( e L φ ′) + A3e L w′ − I 2e L φ

+ A2 − A3e L φ = 0 v

− A1 ∂y ∂y

β

β

y

y

(15)

β

y ∂ ( e L w′) ∂(e L φ )

= 0 A3 − A3 − I 0e L w ∂y ∂y

where a prime and an over-dot represent differentiation with respect to space y and time t respectively. The natural boundary conditions can also be obtained from Eq.(14). The bending moment M , shear force S and axial force F of the BFGT can be expressed as: βy

F = − A0e

L

βy

S = − A3e

L

βy

v′ + A1e L φ ′; M = A1e

βy L

βy

v′ − A2e L φ ′

βy

(16)

w′ + A3e φ L

Assuming harmonic oscillation so that V0 ( y, t ) = V ( y )eiωt ,W0 ( y, t ) = W ( y)eiωt , Φ 0 ( y, t ) = Φ ( y )eiωt

(17)

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

where V ( y ), W ( y ) and Φ( y ) are the amplitudes of V0 ,W0 , Φ 0 , respectively, and ω is the circular frequency. The state space equation of the BFGT can be obtained by substituting Eqs.(12), (13) and (17) into Eq.(15) and Eq.(16), ∂v = Av ∂y

(18)

where v is the state space variable and v = [W , Φ,V , S , M , F ]T . Where M , S , F are the amplitudes of the bending moment, shear force and axial force, respectively. Matrix A can be expressed as: β L  0   A=0  a  5 0 0 

1

β L

0

a1

0

0

0

a2

0

a3

β

0

0

L 0

a6 a7

a7 a5

0

0

−1 0

0 0

 0  a3    a4   0  0 0 

(19)

where a1 = − A3−1 , a2 =

A0 A , a3 = 2 1 A12 − A0 ⋅ A2 A1 − A0 ⋅ A2

(20)

A a4 = 2 2 , a5 = ω 2 I 0 , a6 = ω 2 I 2 , a7 = −ω 2 I1 A1 − A0 ⋅ A2

The special case of a beam with α = 0, β = 0 corresponds to a uniform Timoshenko beam made of homogeneous materials. According to Eq. (18) and using the precise integration method [42], X = eAy ⋅ X0

（21）

where X 0 is the initial state vector of the structure. The six eigenvalues of the matrix A are λi (i = 1 6) . Thus, the analytical solution of the motion differential equation Eq. (15) is: 6

6

6

i =1

i =1

i =1

λi y V ( y ) = ∑ Pe , W ( y ) = ∑ Qi eλi y , Φ ( y ) = ∑ Ri eλi y i

（22）

The following relationships can be derived: Pj = α j R j , Q j = β j R j ( j = 16)

(23)

where αj =

ω 2 I1L + A1λ 2j L − A1βλ j ω 2 I 0 L + A0λ 2j L − A0 βλ j

,β j =

A3 Lλ j

A3 Lλ j2 − A3r β + ω 2 I 0 L

(24)

Similarly the amplitudes of the bending moment M , shear force S and the axial force F can be obtained as: β β   F =  − A0 (V ′ − V ) + A1 (Φ ′ − Φ )  L L   6

（25）

β λ y β = ∑  A0α j − A0λ jα j + A1λ j − A1  ⋅ R j e j L L  j =1 

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

β β   M =  A1 (V ′ − V ) − A2 (Φ′ − Φ )  L L   6 β β λy  = ∑  A1λ jα j − A1α j − A2λ j + A2  R j e j L L  j =1

S = − A3 (W ′ −

β L

（26）

W ) + A3Φ

（27）

6

β λ y   = ∑  − A3β j λ j + A3β j ⋅ + A3  R j e j L  j =1 

The vector R is defined as: R = [ R1 , R2 , R3 , R4 , R5 , R6 ]T

(28)

where T denotes a vector transpose. 2.2 Formation of the dynamic stiffness matrix of the BFGT The dynamic stiffness matrix of the bi-directional functionally graded beam can be formed by the traditional method. The length of the BFGT element is L . The boundary conditions of force and displacement at y = 0 and y = L can be expressed as: At y = 0 : V = Vl ,W = Wl , Φ = Φ l , F = Fl , S = Sl , M = M l (29) At y = L : V = Vr ,W = Wr , Φ = Φ r , F = − Fr , S = − S r , M = − M r (30) The force vector F and the displacement vector d can be expressed as: F = [ Fl , Sl , M l , Fr , Sr , M r ]T , d = [Vl ,Wl , Φ l ,Vr ,Wr , Φ r ]T

(31)

where T denotes a vector transpose. According to Eq. (22), the relationship between R and d is: d = BR (32) where α2 α3 α4 α5 α6   α1  β β2 β3 β4 β5 β 6   1  1 1 1 1 1 B =  λ1L λ2 L λ3 L λ4 L α2 e α 3e α 4e α 5e λ5 L α1e  λ1 L β 2 eλ2 L β 3e λ3 L β 4 eλ4 L β 5e λ5 L  β1e L λ  e1 e λ2 L e λ3 L e λ4 L e λ5 L 

1

α6e



λ6 L 

 λ6 L  β6e  e λ6 L 

According to Eqs. (25)-(27), the relationship between force vector F and vector R is: F = CR Each element of the matrix C can be expressed as: C1 j = C3 j =

β L

β L

(33)

(34)

( A0α j − A1 ) + ( A1 − A0α j ) λ j , C2 j = A3β j  βL − λ j  + A3 ( A2 − A1α j ) + ( A1α j − A2 ) λ j

β  λL C4 j = −  A0α j − A1 + A1 − A0α j λ j  ⋅ e j L 

(

) (

)

（35）

  λL β  C5 j = −  A3 β j  − λ j  + A3  ⋅ e j L    β  λL C6 j = −  A2 − A1α j + A1α j − A2 λ j  ⋅ e j ( j = 1 6) L 

(

) (

)

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

According to Eqs. (32) and (34), the relationship between vector F and vector d can be expressed as: (36) F = K ⋅d where K is the element dynamic stiffness matrix and can be expressed as: K = CB −1

(37) The explicit expressions of the dynamic stiffness matrix of the BFGT can be obtained from Eq. (37). However, it is complicated to compute the dynamic stiffness matrix using Eq. (37). To overcome this shortcoming, a simpler method is proposed to obtain the dynamic stiffness matrix of the beam. The transfer matrix T of the BFGT can be obtained using Eq. (21) as: T = e A ⋅L

(38) According to Eq. (38), the biggest drawback of the transfer matrix method is that the elements of the transfer matrix always contain positive exponential terms. Therefore, the matrix is ill conditioned in the process of computing the natural frequency. However, the dynamic stiffness method can overcome the disadvantages of the transfer matrix method, and because the elements in the dynamic stiffness matrix usually contain only negative exponential terms, there is no overflow or ill conditioned matrix. As the dynamic stiffness method has good computational stability, it is necessary to transform the transfer matrix into the dynamic stiffness matrix. The transfer matrix of the BFGT is  X R   T11 −F  = T  R   21

T12   X L  T22   FL 

(39)

where X R , FR are the state vectors of the displacement and force at the right end of the element and X L , FL are the vectors at the left end of the element. After transformation, Eq.(40) can be obtained from Eq.(39):  FL   K11 F  =  K  R   21

K12   X L  K 22   X R 

(40)

where K11 = − T12−1 ⋅ T11 , K12 = T12−1 , K 21 = − T21 + T22 T12−1T11 K 22 = − T22 T12−1

(41)

Therefore, the element dynamic stiffness matrix K can be obtained using Eq. (41). It is obvious that the dynamic stiffness matrix depends on the frequency. Thus, the dynamic stiffness matrix can be obtained directly from the state space equation through Eqs.(18) and (19). This method is easy to operate and is more general than the traditional method. The method of forming global dynamic stiffness matrix, which can be obtained by assembling element dynamic stiffness matrices, is similar to the finite element method. 2.3Application of the Wittrick-William algorithm For free vibration, according to the boundary conditions, the global dynamic stiffness matrix K (ω ) needs to satisfy: K (ω )d = 0

（42）

where d is the modal shape vector. To obtain the natural frequency of the structure, the dynamic stiffness matrix K (ω ) should satisfy: K (ω ) = 0

(43)

The Wittrick-William algorithm (W-W) [37-38] is an efficient algorithm for solving Eq. (43). This *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

algorithm does not directly compute the natural frequency of the structure. In fact, this method is a counting algorithm. According to the Wittrick-William algorithm, the number of frequencies which are lower than the given value ω* can be given by

{

}

J (ω* ) = J 0 (ω* ) + s K (ω* )

(44)

where J 0 (ω* ) is the total number of frequencies of a individual element with fixed boundary

{

}

conditions that are lower than the given value ω* and s K (ω* ) is the number of negative elements on the leading diagonal of K ∆ . K ∆ is the upper triangular matrix, which can be obtained through

{

}

the method of Gaussian elimination. J , J 0 and s K (ω* ) are integers. The number of natural frequencies in any frequency range can be obtained by the Wittrick-William algorithm. The dichotomy is the simplest way to solve the structural frequencies. For example, to obtain the first order natural frequency, upper bound ωu and lower bound ωl can be obtained by adjusting the test values so that J (ωu ) ≥ i, J (ωl ) ≤ i − 1 (45) Thus, ωl ≤ ω ≤ ωu . In other words, there is only one root in the interval [ωl , ωu ] . A frequency

interval can also be repeatedly divided into two subintervals until the results reach the required accuracy. Because of the low computational efficiency of this method, the author proposes to use a new non-iterative algorithm for solving nonlinear equations, which has been proposed by Yun [43]. This algorithm is suitable for solving nonlinear equations with a single root in the interval (α , β ) . According to literature [43], p* , which is the single root of the nonlinear equation f ( x ) in the interval (α , β ) , can be expressed as: p* =

{

}

β 1 α + β + sgn( f ( x )) ∫ sgn( f ( x ))dx α 2

(46)

According to Eq. (46), the following algorithm can be obtained: 1. p* is the single root of nonlinear equation f ( x ) in the interval ( a, b) . First, the interval is divided into N segments with δ =

(b − a ) . Using numerical integration, p* can be expressed as: 2N N −1

p∗ = p + δ ⋅ F (a ) ⋅ ∑ F ( p + (2 j − N )δ )

(47)

j =1

where p=

a+b , F ( x ) = sgn ( f ( x ) ) 2

(48)

2. Newton’s method takes p* as the initial value to compute the root of f ( x ) . The convergence of the method is demonstrated in detail in the literature [43]. 2.4 Frequency response function of the BFGT The frequency response function matrix H(ω ) of the BFGT is the inverse of the global *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

dynamic stiffness matrix D(ω ) : [H(ω )] = [D(ω )]−1

(49)

2.5 Dynamic response of the BFGT under a harmonic moving load The harmonic moving load can be expressed as: P(t ) = P0 cos(Ωt )δ ( x − v ⋅ t ) (50) where P0 is the amplitude of the harmonic load, δ () is the Dirac function and Ω is excitation frequency of the moving load. Thus, the potential energy of the moving load can be written as: L

U ext = − ∫ P(t ) ⋅ w( x , t )dx

(51)

0

According to the mode superposition principle, the displacement response w, v, φ can be expressed as: N

w( y , t ) = ∑ wi ( x ) Ri (t ) i =1 N

v( x , t ) = ∑ vi ( x ) Si (t )

(52)

i =1 N

φ ( x, t ) = ∑ φi ( x )Ti (t ) i =1

where Ri (t ) Si (t ) Ti (t ) are unknown time-dependent coefficients and wi ( x ), vi ( x ), φi ( x ) are mode shapes that need to satisfy the geometric boundary conditions. wi ( x ), vi ( x ), φi ( x ) can be obtained from Eq.(22). The motion equations of the BFGT can be derived using the Lagrange equations as follows: d  ∂T  dt  ∂R
i

 ∂U ∂U ext + = 0(i = 1 N ) + ∂Ri  ∂Ri

(53)

d  ∂T  ∂U ∂U ext + = 0(i = 1 N )  + dt  ∂S
i  ∂Si ∂Si d  ∂T  ∂U ∂U ext + = 0(i = 1 N )  + dt  ∂T
i  ∂Ti ∂Ti

Thus,

( t )} [K ] [0] [K ]  {R (t )} {F( t )} [0 ]  {R [M11 ] [0] 11 13       [0] [ M ] [ M ]  {

S(t )}  +  [0] [ K 22 ] [ K 23 ]  {S(t )}  =  {0}  22 23     

(t )} [ K ] [ K ] [ K ] {T( t )}   {0}   [0] [M 32 ] [ M 33 ] {T    32 33     31

(54)

where βy

L

[K11 ]ij = ∫ A3e

βy βy βy L ∂wi ∂w j ∂wi L L ∂v ∂v ∂v ∂φ dy ,[K13 ]ij = − ∫ A3e L ⋅ ⋅ φ j dy [K 22 ]ij = ∫ A0e L i ⋅ j dy ,[ K 23 ]ij = − ∫ A1e L i ⋅ j dy 0 0 0 ∂y ∂y ∂y ∂y ∂y ∂y ∂y

L

0

L

βy

[K 31 ]ij = − ∫ A3e

L

0

φi

∂w j ∂y

dy

βy βy L ∂φi ∂φ j  [ Κ 33 ]ij = ∫  A3e L φi ⋅ φ j + A2e L ⋅ dy 0  ∂y ∂y  

βy

L

[M11 ]ij = ∫ I 0 e 0

L

L

L

wi w j dy [M 22 ]ij = ∫ I 0 e 0

βy

[M 32 ]ij = − ∫ I1e L φi v j dy {F( t )}i = 0

βy L

L

vi v j dy [ M 23 ]ij = − ∫ I1e 0

L

βy L

L

[K 32 ]ij = − ∫ A1e 0

L

βy L

∂φi ∂v j dy ∂y ∂y

βy

viφ j dy [M33 ]ij = ∫ I 2 e L φiφ j dy 0

∫0 P(t )wi dt = P0 cos(Ω ⋅ vt ) ⋅ wi (vt ) (i, j = 1, 2, 3 N )

（55）

Eq. (54) can also be written as:

+ Kx = F Mx

(56)

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

where T

T

(57)

x = {{R (t )},{S(t )},{T(t )}} F = {{F(t )},{0},{0}}

Eq. (56) can also be written as: u
= H ⋅ u + f

(58)

 [0] M −1  0 u = {x, Mx
}T , H =  ,f =    F  − K [0] 

(59)

where

According to the precise integration method [42], thus, t

u( t ) = e H ( t −t0 ) u( t0 ) + ∫ eH ( t −τ ) f (τ ) dτ

(60)

t0

The discrete form of the Eq. (60) is u(( j + 1) h ) = eHh u( jh ) + ∫

( j +1) h H (( j +1) h −τ ) jh

e

(61)

f (τ ) dτ

Where h is integral time step. f (τ ) is the constant in the time interval [ jh ,( j + 1)h ] ,thus, u(( j + 1)h) = eHh u( jh) + H−1 ( eHh − I)f ( jh )

(62)

Where eHh can be computed as described in literature [42]. Therefore, the dynamic response of the BFGT under a harmonic moving load can be computed using Eqs. (62) and (52).

3. Numerical results 3.1 Validation Comparison is needed to verify the correctness of the method. The results obtained by this method must be compared with the previous literature. The derivation in this paper is suitable for any form of distribution along the thickness direction including the power law [41] distribution. To verify the results of this method, two cases are discussed: First case: the geometrical and material properties are as follows: L = 0.4m , Elb = 163.8Gpa Glb = 78.8Gpa , ρlb = 7850kg / m3 , shear coefficient k s = 5 / 6 , section area A = 1.6 × 10−3 m 2 , section

moment of inertia I = 2.56 × 10−6 m 4 ,and gradient parameter α = 0 . Thus, the material properties only vary along the axial direction. The natural frequencies under different boundary conditions are compared with Ref[19], as shown in Table 1. The non-dimensional frequency, as defined in Ref[19], is:

λ = ω L2

ρ lb A

(63)

Elb I

Ref[19] obtained analytical solutions of exponentially non-uniform functionally graded Timoshenko beams, and it is observed that the present results are in good agreement with analytical values in table 1. Table 1 Comparison of dimensionless natural frequencies with axially exponential gradient (α = 0) β 0

Mode

C −C

number

Ref[19]

1

13.8348

S−S

S−S

C−F

C−F

Present

Ref[19]

Present

Ref[19]

Present

13.834758

8.3874

8.387358

0.6257

0.625733

C −C

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

1

2

2

28.5179

28.517925

25.3460

25.345880

16.7920

16.791956

3

45.6660

45.665951

44.1266

44.126572

33.8149

33.814870

4

61.8621

61.862051

56.6139

56.613648

51.5214

51.521414

1

13.9007

13.900718

8.2533

8.253299

0.6273

0.627255

2

28.6185

28.618458

25.3489

25.348854

16.8358

16.835847

3

45.7468

45.746811

44.1083

44.108334

33.8702

33.870186

4

61.9873

61.987329

56.8471

56.847134

51.5338

51.533772

1

14.1133

14.113294

7.8644

7.864356

0.6265

0.626454

2

28.9175

28.917533

25.3649

25.364944

16.9805

16.980545

3

45.9881

45.988096

44.0686

44.068596

34.0383

34.038303

4

62.3482

62.348191

57.5072

57.507238

51.5878

51.587773

Second case: The derivation in this paper is suitable for any form of distribution in the thickness direction. Ref[41] obtained the exact natural frequencies of functionally graded Timoshenko beams. In Ref[41], the material properties vary continuously in the thickness direction according to the power law distribution, but the material properties do not vary along the axial direction. The geometrical and material properties of the beam are as follows: Elb = 210Gpa , ρlb = 7800kg / m3 , Elt = 390Gpa , ρ lt = 3960kg / m3 , ν = 0.3 , h = 0.1m , b = 0.1m . k is the power law index defined

in Ref[57]. The shear correction factor has the value of ks = 5 / 6 and non-dimensional frequency is defined as: λ=

ω L2 ( ρlb / Elb )

(64)

h

By comparing the calculated frequency values with the known results [41], it is observed that present dimensionless frequency coincide exactly with the previous results according to table2. Table 2

Comparsion of dimensionless natural frequencies ( β = 0) under SS boundary condition k = 0.1

k = 0.1

k = 0.5

k =1

k =2

L/h

method

5

Ref[41]

4.784

4.529

4.059

3.689

3.390

present

4.784023

4.52912

4.059011

3.68891

3.39045

10

Ref[41]

16.652

15.770

14.128

12.818

11.740

present

16.65169

15.77037

14.12844

12.81768

11.73961

20

Ref[41]

28.189

26.780

24.022

21.621

19.479

present

28.18882

26.78042

24.02192

21.62066

19.47922

30

Ref[41]

31.924

30.239

27.085

24.531

22.387

present

31.92443

30.23861

27.08543

24.53057

22.38664

3.2 Dynamic characteristics analysis of the BFGT The geometrical and material properties of a bi-directional functionally graded Timoshenko beam of unit width are as follows: Elb = 71Gpa , ρlb = 2780kg / m3 , ν = 0.3 , h = 0.05m , L = 5m and shear coefficient ks = 5 / 6 . Non-dimensional frequency is defined as:

λ=

ω L2 ( ρlb / Elb ) h

(65)

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

Exact solutions of the first-order natural frequency of the BFGT under different boundary conditions are shown in table 3-6, Where S , C and F denote simple, clamped and free ends, respectively. Table 3 Dimensionless fundamental frequency with CC boundary conditions λ

α =0

α =2

α =4

α =6

α =8

α = 10

β =0

6.4541

5.8729

4.6643

3.5570

2.7661

2.2321

β =2

6.6168

6.0210

4.7820

3.6467

2.8359

2.2884

β =4

7.1506

6.5068

5.1679

3.9411

3.0648

2.4731

β =6

8.1620

7.4273

5.8990

4.4987

3.4985

2.8231

β =8

9.7532

8.8753

7.0493

5.3761

4.1808

3.3737

β = 10

11.9796

10.9015

8.6589

6.6037

5.1356

4.1442

Table 4 Dimensionless fundamental frequency with SS boundary conditions λ

α =0

α =2

α =4

α =6

α =8

α = 10

β =0

2.8486

2.9203

2.9251

2.6919

2.3408

2.0075

β =2

2.7382

2.7853

2.7303

2.4480

2.0859

1.7659

β =4

2.4272

2.4310

2.2927

1.9755

1.6387

1.3662

β =6

1.9746

1.9525

1.7894

1.5031

1.2280

1.0154

β =8

1.4675

1.4401

1.2994

1.0777

0.8741

0.7201

β = 10

0.9960

0.9740

0.8727

0.7198

0.5820

0.4787

Table 5 Dimensionless fundamental frequency with CS boundary conditions λ

α =0

α =2

α =4

α =6

α =8

α = 10

β =0

4.4491

4.1417

3.5071

2.8992

2.4096

2.0322

β =2

5.1133

4.7703

4.0275

3.2681

2.6499

2.1905

β =4

5.9953

5.5753

4.6472

3.6989

2.9516

2.4152

β =6

7.2641

6.7180

5.5203

4.3295

3.4213

2.7840

β =8

9.0944

8.3267

6.7639

5.2512

4.1245

3.3454

β = 10

11.4217

10.4697

8.4377

6.5081

5.0928

4.1227

Table 6 Dimensionless fundamental frequency with CF boundary conditions λ

α =0

α =2

α =4

α =6

α =8

β =0

1.0149

0.9234

0.7332

0.5591

0.4348

0.3508

β =2

0.5313

0.4834

0.3838

0.2927

0.2276

0.1836

β =4

0.2614

0.2378

0.1888

0.1440

0.1119

0.0903

β =6

0.1222

0.1112

0.0883

0.0673

0.0523

0.0422

β =8

0.0549

0.0500

0.0397

0.0302

0.0235

0.0190

β = 10

0.0239

0.0218

0.0173

0.0132

0.0103

0.0083

α = 10

The effect of gradient parameters α , β on the first-order natural frequency of the BFGT is *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

compared in Table 3-6. For CC boundary conditions and CS boundary conditions, increasing the axial gradient parameter β increases the fundamental frequency of the BFGT. However, increasing the gradient parameter α along the thickness direction makes the structure's fundamental frequency decrease. In contrast, for CF boundary conditions, the fundamental frequency values decrease as the material gradient parameters α , β increase. For SS boundary conditions, with the increase of the axial gradient parameter β , the fundamental frequency of the BFGT decreases, but the influence of gradient parameter α on the fundamental frequency is more complex and can be affected by the axial gradient parameter β . As observed from the tabulated results, the fundamental frequency of the BFGT is significantly influenced by the axial gradient parameter β . (a)C-S

(b) C-F

350

200 α=0, β =0 α=0, β =5

250

α =0, β =0 α =0, β =5 α =0, β =-5 α =5, β =0 α =10,β =0

180

α=0, β =-5 α=5, β =0 α=10, β =0

160 circular frequency(rad/s)

circular frequency(rad/s)

300

200 150

100

140 120 100 80 60 40

50 20

0 20

30

40

50 60 70 slenderness ratio L/h

80

90

0 20

100

(c)C-C

40

50 60 70 slenderness ratio L/h

80

90

100

(d) S-S 400

160

300 250 200 150 100

α =0, β =-5 α =5, β =0 α =10,β =0

120 100 80

60

40

50 0 20

α =0, β =0 α =0, β =5

140

circular frequency(rad/s)

α =0, β =0 α =0, β =5 α =0, β =-5 α =5, β =0 α =10, β =0

350

circular frequency(rad/s)

30

30

40

50 60 70 slenderness ratio L/h

80

90

100

20 20

30

40

50 60 70 slenderness ratio L/h

80

90

100

Fig2 Fundamental frequency of BFGT beams with different slenderness ratios: (a) clamped-simple (b) clamped-free(c) clamped-clamped (d)simple-simple

Fig. 2 shows the influence of the slenderness ratio L / h on the fundamental frequency under different boundary conditions when the length of the BFGT is constant, where the value of h varies and the other material parameters are the same as before. Fig. 2 shows that the fundamental frequency of the BFGT increases as the slenderness ratio decreases. When the slenderness ratio is less than 50, the shear effect will have a significant impact on the natural frequency of the BFGT. For the CC and SS boundary conditions. The curves corresponding to β = 5 and β = −5 are completely coincident due to the symmetric boundary conditions. For the CS and CF boundarys, there is a significant difference between the curves corresponding to to β = 5 and β = −5 because of the asymmetric boundary conditions. Fig. 3 illustrates the influence of gradient parameters α , β on the normalized mode shapes of the BFGT. Fig.a shows the influence of the gradient parameters on the first order mode and the third order mode under the SS boundary condition. Fig.b and Fig.c *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

show the influence of gradient parameters on the first order mode and the second order mode under the CC and CF boundary conditions. It is obvious that gradient parameter β has a significant influence on mode shapes compared to gradient parameter α , and the variation of gradient parameter β can change the location of the node of the mode shape. The mode shapes are symmetrical when β = 0 , while the mode shapes are asymmetric when β ≠ 0 according to Fig. 3. (a)

(b)

(c)

1.2

1 α =0,β =0 α =5,β =0

α =0,β =0 α =5,β =0

1

α =0,β =10

0.8

0.6

0.4

0.2

0

0.6

0.4 0.353 0.2

0.5

1

1.5

2

2.5 3 location(m)

3.5

4

4.5

-0.2 0

5

normalized mode shapes

-0.5 α =0,β =0 α =5,β =0 α =10,β =0 α =0,β =5

-1

2.5212 2.5214 2.5216

0.5 0.4 α =0,β =0

0.3 0.2

0.3515

0.1

0.936 0.937 0.938 0.5

1

1.5

2

2.5 3 location(m)

3.5

4

4.5

0

5

α =5,β =0

0

α =10,β =0 α =0,β =5 α =0,β =5

0.5

1

1.5

2

2.5 3 location(m)

0.5

1

1.5

2

2.5 3 location(m)

3.5

4

4.5

4

4.5

5

0.5

0

-0.6537 -0.5 -0.6537 -0.6537 -0.6537 2.9705 2.971 2.9715 2.972 2.9725 -1 location(m)

0

-0.5

0.5052 0.5052 0.5051 0.5051 0.5051 1.133 1.134 1.135

α =0,β =0 α =5,β =0 α =10,β =0 α =0,β =5

-1

α =0,β =10

-1.5 0

3.5

1 α =0,β =0 α =5,β =0 α =10,β =0 α =0,β =5 α =0,β =10

0.5

0

0.6

0.352

1

0.5

0.7

0.3446 0.3446 0.3446 0.3446 0.3445

0.3525

0

1

normalized mode shapes

0.8

0.8

normalized mode shapes

-0.2 0

0.9

α =10,β =0 α =0,β =5 α =0,β =10

α =10,β =0 α =0,β =5

normalized mode shapes

normalized mode shapes

1

normalized mode shapes

1.2

α =0,β =10

5

-1.5 0

0.5

1

1.5

2

2.5 3 location(m)

3.5

4

4.5

5

-1.5 0

0.5

1

1.5

2

2.5 3 location(m)

3.5

4

4.5

5

Fig3 Influence of gradient parameters α , β on BFGT mode shapes under different boundary conditions (a) simple-simple (b)clamped-clamped (c) clamped-free

Fig. 4 depicts the influence of gradient parameters α , β on the transverse displacement origin FRF at the midspan ( x = L / 2) of the BFGT under different boundary conditions. The logarithmic amplitude frequency diagram of the FRF can not only reflect amplitude-frequency characteristics of the system in a wide frequency band, but can also represent the prominent features of te amplitude-frequency curve in both the low and high bands of frequency. According to Fig.a, within the band 0-150Hz, the curve in the low frequency band is approximately a straight horizontal line. Because the response curve in the low frequency band depends on the stiffness of the system, it is also known as the stiffness line, and the curve in the high frequency band approximates an inclined line that primarily depends on the mass of the system and is known as mass line. Thus, for the CC and SS boundary conditions, increasing the gradient parameters α , β can increase the mass and stiffness of the BFGT. For the CF boundary condition, as gradient parameter β increases, the stiffness of the structure decreases, but the mass increases. However, increasing the gradient parameter α makes both the stiffness and mass of the BFGT increase. For the CC and SS boundary conditions, curve mutation of the FRF near the second-order frequency can be caused by variation of the gradient parameter β . This is because the midspan of the BFGT is the modal node when β = 0 , while the node position of the second mode is changed due to β ≠ 0 . Fig.4 shows that the frequency difference between the adjacent resonance and anti-resonance frequencies can be significantly affected by the gradient parameter β . (a)

(b)

(c)

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

-60

-60

-80

-60

-80

-80

-100

-100

-120

-120

-90 -100

-100

-110

-160

-140

α =0,β =0 α =0,β =1 α =0,β =2 α =0,β =3

-200 -220

0

50

100

150 200 250 circular frequency(rad/s)

300

350

-180

α =0,β =2

-200

50

100

150 200 250 circular frequency(rad/s)

300

350

α =0,β =3

-220 0

400

-100

-110 -100

-120 0

20

40

-140 -160

-200

50

100

150 200 250 circular frequency(rad/s)

300

350

-240

400

α =2,β =0 α =3,β =0 α =4,β =0

20 40 60 80

-140

-180

α =2,β =0 α =3,β =0

-220

α=4,β =0

-220 0

α =0,β =0 α =1,β =0

-160 α =1,β =0

α=2,β =0 α=3,β =0

400

α =0,β =0

-180

α=1,β =0

350

-120

-120

60 FRF(db)

FRF(db)

-120

α=0,β =0

300

-110

-100

-100

-160

150 200 250 circular frequency(rad/s)

-80

-90

-80

-140

100

-100

-60

-120

50

-60

-40

-200

α =0,β =0 α =0,β =1 α =0,β =2

-200

α =0,β =3

-220 0

400

-80

FRF(db)

-180

α =0,β =0 α =0,β =1

-60

-180

0 20 40 60

-140 -160

-160

-180

-240

FRF(db)

-140

FRF(db)

FRF(db)

-120

-200

α =4,β =0

0

50

100

150 200 250 circular frequency(rad/s)

300

350

-220 0

400

50

100

150 200 250 circular frequency(rad/s)

300

350

400

Fig4 Frequency response function of the BFGT (a) clamped-clamped (b) clamped-free (c) simple-simple

Fig. 5 shows the effect of the speed of a moving harmonic load on the maximum transverse displacement, velocity and acceleration at midspan ( x = L / 2) under the SS boundary condition. The magnitude of the moving load is 10N . When the excitation frequency Ω = 0rad / s , it is clear that there exists a critical velocity v = vc that maximizes the deflection. Because increasing the gradient parameters α , β can increase the stiffness of the BFGT, the maximum deflection at mid span decreases as the gradient parameters α , β increase. In Fig.b, when the speed of moving load exceeds 30m / s , the variation of the maximum speed at mid span is slow and approximates a constant, and this constant decreases with the increase of the parameters α , β . (a) (b) (c) -5

-3

0.9

5.5

0.8

5 4.5 4 3.5

α =0,β =0 α =0,β =1

3

2

0

10

20 30 Load velocity,v(m/s)

40

0.7 0.6 0.5 0.4 α =0,β=0

0.3 0.2

α =0,β=2 α =0,β=3 α =0,β=4

0.1 0

50

-5

0

x 10

1

x 10

α =0,β =0

5

10

15

20 25 30 Load velocity,v(m/s)

35

40

45

4

3

2

α =0,β =0

0.01

0

50

α =0,β =1 α =0,β =2 α =0,β =3 α =0,β =4

0

5

10

15

20 25 30 Load velocity,v(m/s)

35

40

45

50

0.04 α =0,β =0 α =1,β =0

0.035

0.8

α =0,β =0 α =1,β =0

0.7

α =2,β =0 α =3,β =0 α =4,β =0

0.6

α =2,β =0 2 Maximum Acceleration(m/s )

Maximum velocity(m/s)

α =4,β =0

5

0.02 0.015

-3

0.9

α =2,β =0 α =3,β =0

0.03 0.025

0.005

α =1,β =0

6 Maximum displacement(m)

0.04 0.035

α =0,β=1

α =0,β =2 α =0,β =3 α =0,β =4

2.5

7

x 10

1

6

Maximum Acceleration(m/s2)

x 10

Max imum velocity(m/s)

Maximum displacement(m)

6.5

0.5 0.4 0.3

0.03

α =3,β =0 α =4,β =0

0.025 0.02 0.015 0.01

0.2 1

0

0.005

0.1 0

5

10

15

20 25 30 Load velocity ,v(m/s)

35

40

45

50

0

0

5

10

15

20 25 30 Load velocity,v(m/s)

35

40

45

50

0

0

5

10

15

20 25 30 Load velocity,v(m/s)

35

40

45

50

Fig5 Variation of the maximum dynamic response with the moving load velocity under SS boundary conditions: (a) maximum displacement (b) maximum velocity (c) maximum acceleration

Fig. 6 is the displacement dynamic response at the midspan of the bi-directional functionally graded Timoshenko beam. The influences of the moving velocity, excitation frequency and gradient parameters on the response of the BFGT are compared in Fig. 6, where the vertical coordinate is the normalized deformation. When the excitation frequency is close to the *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

fundamental frequency of the BFGT, the amplitude of the structure increases significantly, which is called the resonance phenomenon. Since the variation of parameters α , β will directly affect the fundamental frequency of the BFGT, it is obvious that the parameters α , β have a great impact on the dynamic response. (a) v = 1m / s, Ω = 0rad / s

(b) v = 1m / s, Ω = 20rad / s

1.2

0.5

0.6

0.4 0.8 0.2

0.7

1

0.5

1

0

-0.5 0.4 α =0,β =0

-1

0.6 1.5

1.5

2

3

3.5

4

4.5

-1.5 0

5

(d) v = 5m / s, Ω = 0rad / s

0.5

1

1.5

0

-0.5

-1

0.2

α =10,β =0 α =0,β =10

2

2.5 time(s)

Normalized displacement

0.8

Normalized displacement

0 4.2 2

2.5 time(s)

3

3.5

4.4 4

4.5

-1.5 0

5

(e) v = 5m / s, Ω = 20rad / s α =0,β =0 α =10,β =0 α =0,β =10

1

1.5

2

2.5 time(s)

3

3.5

4

4.5

5

1 0.8

0.5 Normalized displacement

1

0.5

(f) v = 5m / s, Ω = 40rad / s

1

1.2

α =0,β =0 α =10,β =0 α =0,β =10

4.6

0.8

0.6

0.4

0.6 Normalized displacement

Normalized displacement

0.5

α =0,β =10

0

Normalized displacement

1

α =0,β =0 α =10,β =0

1

-0.2 0

(c) v = 1m / s, Ω = 40rad / s

1

0

-0.5

0.4 0.2 0 -0.2 -0.4

0.2

-1

-0.6

α =0,β =0 α =10,β =0

0

α =0,β =0 α =10,β =0

-0.8

α =0,β =10

-0.2 0

-1.5 0 0.1

0.2

0.3

0.4

0.5 time(s)

0.6

0.7

0.8

0.9

0.1

0.2

0.3

0.4

1

(g) v = 10m / s, Ω = 0rad / s

0.5 time(s)

0.6

0.7

0.8

0.9

α =0,β =10

-1 0

1

(h) v = 10m / s, Ω = 20rad / s

1

0.1

0.2

0.3

0.4

0.5 time(s)

0.6

0.7

0.8

0.9

1

(i) v = 10m / s, Ω = 40rad / s

1

1

0.5

0.5

0.4

0.2

0 α=0,β =0 α=10,β =0 α=0,β =10

-0.2

-0.4 0

0.05

0.1

0.15

0.2

0.25 time(s)

0.3

0.35

0.4

0.45

Normalized displacement

0.6

Normalized displacement

Normalized displacement

0.8

0

-0.5

-1

0

-0.5

-1

α =0,β =0

α =0,β =0 α =10,β=0

α =10,β =0 α =0,β =10

0.5

-1.5 0

0.05

0.1

0.15

0.2

0.25 time(s)

0.3

0.35

0.4

0.45

α =0,β =10

0.5

-1.5 0

0.05

0.1

0.15

0.2

0.25 time(s)

0.3

0.35

0.4

0.45

0.5

Fig6 Time history of the normalized dynamic deflection

4. Conclusion In this paper, the differential equations of a bi-functional graded Timoshenko beam are established using the Hamilton’s principle. The exact dynamic stiffness matrix of the BFGT is formed and the frequency is solved by the Wittrick-William algorithm and a non-iterative algorithm. 1. In this paper, the influence of gradient parameters α , β on the fundamental frequency of the BFGT under different boundary conditions is analysed. Increasing gradient parameter α can decrease the fundamental frequency of the BFGT, while the influence of axial gradient parameter β on the fundamental frequency of the BFGT is related to the boundary conditions. 2. The slenderness ratio has a great impact on the fundamental frequency of the BFGT. When the structural slenderness ratio decreases, the shear effect of the structure increases significantly 3. This paper analyses the influence of gradient parameters α , β on the frequency response function of the BFGT, and the stiffness and mass characteristics of the BFGT are analysed using the curves of the FRF under different boundary conditions. *Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

4. The influence of the gradient parameters α , β on the modal shape of the BFGT is analysed. The results show that the gradient parameter α has little influence on the vibration modes of the structure, while the influence of the axial gradient parameter β on the vibration modes is relatively significant. 5. Differential equations of the BFGT under harmonic moving loads are derived by means of Lagrange equations, and the dynamic displacement response of the BFGT is solved using the modal superposition method. The maximum transverse displacement, velocity and acceleration are greatly influenced by parameters α , β under the SS boundary condition and the parameters have a great impact on the time history of deflection. 6. The method of forming the dynamic stiffness matrix directly from the state space differential equations is versatile, and the derivation in this paper is not only suitable for material properties varied exponentially along the thickness direction, but also for other forms of distribution.

References [1] Ichikawa K,editor.Functionally graded materials in the 21 st century:a workshop on trends and forecasts.Japan:Kluwer Academic Publishers:2000 [2] Yang J,Xiang HJ.Thermo-electro-mechanical characteristics of functionally graded piezoelectric actuators.Smart Mater Struct 2007;16:784-97 [3] BirmanV,Byrd LW.Modeling and analysis of functionally graded materials and structures.Appl Mech Rev 2007;60:195-216 [4] Analytical solution of a cantilever functionally graded beam.Compos Sci Technol 2007;61:689-696 [5] Li X-F.A unified approach for analyzing static and dynamic behaviors of functionally graded Timoshenko and Euler-Bernoulli beams.J Sound Vib 2008;318:1210-29 [6] Li X-F A higher-order theory for static and dynamic analyses of functionally graded beams.Arch Appl Mech 2010;80: 1197 -1212 [7] Sankar BV.An elasticicity solution for functionally graded beam.Compos Sci Technol 2001;61:689-696 [8] Li X-F.A unified approach for analyzing static and dynamic behaviors of functionally graded Timoshenko and Euler-Bernoulli beams.Journal of Sound and vibration.2008;318:1210-1229 [9]Simsek M,Kocaturk T.Free and forced vibration of a functionally graded beam subjected to a concentrated moving harmonic load. Composite Structures.2009;90:465-473 [10]Pradhan KK,Chakraverty S.Free vibration of Euler and Timoshenko functionally graded beams by Rayleigh-Ritz method.Composite-Part B.2013;51:175-184 [11] Elishakoff I.Eigenvalues of inhomogeneous structures: unusual closed-form solutions.Boca Raton:CRC Press;2005 [12] Ece MC,Aydogdu M,TaskinV.Vibration of a variable cross-section beam.Mech Res Commun 2007;34:78-84 [13] Atmane HA, Tounsi A,Meftah SA,Belhadj HA.Free vibration behavior of exponentially functionally graded beams with varying cross-section.J Vib Control 2011;17:311-318 [14] Wang CY,Wang CM.Exact vibration solutions for a class of nonuniform beams.J Eng Mech.2013;139:928-931 [15] Wang CY, Wang CM.Exact vibration solution for exponentially tapered cantilever with tip mass.Vib Acoust.2012;134 [16] Rajasekaran S. Differential transformation and differential quadrature methods for centrifugally stiffened axially functionally graded tapered beams.Int J Mech Sci 2013;74:15-31 [17] Esmailzadeh E,Ohadi A.Vibration and stability analysis of non-uniform Timoshenko beams under axial and distributed tangential load.Journal of Sound and Vibration 2000;236:443-456 [18] Huang Y,Yang L-E,Luo Q-Z.Free vibration of axially functionally graded Timoshenko beams with non-uniform cross-uniform cross-section.Composite Part B Eng 2012;45:1493-1498 [19] A-Y Tang,J-X Wu.Exact frequency equations of free vibration of exponentially non-uniform functionally graded Timoshenko beams.International Journal of Mechanical

Sciences, 2014;89:1-11

[20]M.Simsek,Fundamental frequency analysis of functionally graded beams by using different high-order beam-theories, Nucl,Eng. Des.240(2010):697-705 [21] Salamat-talab M, Nateghi A, Torabi J. Static and dynamic analysis of thirdorder shear deformation FG micro beam based on modified couple stresstheory. Int J Mech Sci 2012;57(1):63–73. [22] Vo TP, Thai HT, Nguyen TK, Maheri A, Lee J. Finite element model for vibration and buckling of functionally graded sandwich beams based on a refined shear deformation theory. Eng Struct 2014;64:12–22. [23]Zhang B, He Y, Liu D, Gan Z, Shen L. Size-dependent functionally graded beam model based on an improved third-order shear deformation theory. Eur J Mech-A/Solids 2014;47:211–30 [24]Simsek M, Reddy JN. A unified higher order beam theory for buckling of a functionally graded microbeam embedded in elastic

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]

medium using modified couple stress theory. Compos Struct 2013;101:47–58 [25] Ke LL, Wang YS. Size effect on dynamic stability of functionally graded microbeams based on a modified couple stress theory. Compos Struct 2011;93(2):342–50. [26] Nemat-Alla M. Reduction of thermal stresses by developing two-dimensional functionally graded materials. Int J Solids Struct 2003;40:7339–56. [27]Lü CF, Chen WQ, Xu RQ, Lim CW. Semi-analytical elasticity solutions for bidirectional functionally graded beams. Int J Solids Struct 2008;45:258–75. [28] Nie GJ, Zhong Z. Semi-analytical solution for three dimensional vibration of functionally graded circular plates. Computer Methods in Applied Mechanics and Engineering, 2007;196(49-52): 4901-4910 [29] Nie GJ, Zhong Z. Axisymmetric bending of two directional functionally graded circular and annular plates. Acta Mechanica Solida Sinica, 2007; 20(4):289-295 [30] Nie GJ, Zhong Z. Vibration analysis of functionally graded annular sectorial plates with simply supported radial edges. Composite Structures, 2008; 84(2): 167-176 [31] Nie GJ, Zhong Z. Asymmetric bending of bi-directional functionally graded circular plates. In: Fan J H, Chen H B, eds. Advances in Heterogeneous Material Mechanics 2008, Proceedings of the Second International Conference on Heterogeneous Material Mechanics, Huangshan(China), 2008-06-03-08. Lancaster, Pennsylvania:DEStech Publications Inc, 2008. 1013-1016 [32] Mesut Simsek,Bi-directional functionally graded materials(BDFGMs) for free and forced vibration of Timoshenko beams with various boundary conditions.Composite Structures 133(2015):968-978 [33]Banerjee J R.Coupled Bending torsional Dynamic Stiffness Matrix for Beam Elements.International Journal for Numerical Methods in Engineering 1989;28:1283-1289 [34] Banerjee J R.Fisher S A.Coupled Bending-torsional Dynamic Stiffness Matrix for Axially Loaded Beam Elements.International Journal for Numerical Methods in Engineering.1992;33:739-751. [35] Banerjee J R,Williams F W.Coupled Bending torsional Dynamic Stiffness Matrix for Timoshenko Beam Elements. Computers and Styructure. 1992;42:301-310 [36] Li Jun,Shen Rongying,Hua Hongxing,Coupled Bending and Torsional Vibration of Axially Loaded Bernoulli Euler Beams Including Warping Effects.Applied Acoustics,2004;65:153-170 [37] Banerjee JR,Williams F W.Coupled Bending torsional Dynamic Stiffness Matrix of an Axially Loaded Timoshenko Beam Element.International Journal of Solids and Structures.1994;31:749-762 [38] Banerjee J R,Williams FW.Exact Bernoulli Euler Dynamic Stiffness Matrix for a Range of Tapered Beams.International Journal for Numerical Methods in Engineering.1985;21:2289-2302 [39] Leung A Y T.Zhou W E.Dynamic Stiffness Analysis of Non-uniform Timoshenko Beams Journal of Sound and Vibration.1995;181:447-456 [40]Kolousek V.Anwendung des Gesetzes der virtuellen Ver schiebungen und des Reziprozitatssatzes in der Stabwerksdy nam ik Ingenieur Archiv.1941;21:363-340 [41]H.Su,J.R.Banerjee.Development of dynamic stiffness method for free vibration of functionally graded Timoshenko beams. Comput Struct 2015;147:107-116 [42]X.Q.Zhu,S.S.Law.Precise time-step integeration for the dynamic response of continuous beam under moving loads.Journal of sound and vibration 2001;240(5):962-970 [43] Beong ln Yun,Iterative methods for solving nonlinear equations with finitely many roots in an interval.Journal of Computational and Applied Mathematics.2012,236:3308-3318

*Corresponding author at: School of Aeronautic Science and Engineering , Beihang University, Beijing 100191, China . E-mail address:[email protected]