Influences of hygrothermal environment and installation mode on vibration characteristics of a rotating laminated composite beam

Mechanical Systems and Signal Processing 91 (2017) 23–40

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Influences of hygrothermal environment and installation mode on vibration characteristics of a rotating laminated composite beam Y. Qin, Y.H. Li ⇑ School of Mechanics and Engineering, Southwest Jiaotong University, Chengdu 610031, PR China

a r t i c l e

i n f o

Article history: Received 15 September 2016 Received in revised form 16 December 2016 Accepted 31 December 2016

Keywords: Flapwise free vibration Rotating laminated composite beam Green’s function Hygrothermal environment

a b s t r a c t Flapwise vibration characteristics of a rotating laminated composite beam under hygrothermal environment are studied in this paper. The governing equation of the flapwise vibration for the rotating beam under hygrothermal environment is obtained applying D’Alembert’s principle, in which the design parameters such as setting angle and pitch angle are considered. A numerical method based on Green’s function to handle the vibration characteristics problem of a rotating composite beam is developed. Influences of the hub radius, rotating speed, pitch angle, setting angle, fiber orientation angles, temperature and humidity on natural frequencies of the rotating laminated composite beam are discussed. Results indicate that: (1) the rotating speed, hub radius, setting angle, pitch angle and fiber orientation angle have more prominent effects on natural frequencies than the temperature variation and moisture concentration do; (2) the influences of some parameters such as the hub radius, rotating speed, setting angle and pitch angle on the vibration characteristics are not independent, but coupled with each other. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, composite structures have been used extensively in aerospace and other engineering applications. The composite materials, which provide the structural configuration possibilities, are vital to affect the dynamic behavior of rotating structures. The studies on the composite structures have been attracting attentions by more and more researchers [1–4] in recent years. These structures usually operate in complex environmental conditions, such as hygrothermal environment. On one hand, composite materials will be affected by aging behaviors resulted from hygrothermal environment, which can change the mechanical and physical behaviors of composites. On the other, the variation of temperature and moisture concentration can induce the change of material properties of the composite structures, which naturally reduces the stiffness and strength of the structures. Therefore the vibration analysis of such beam structures has attracted more and more attentions. In order to explore this type of structural systems, a good knowledge of their vibration characteristics is essential. Numerous previous researches relating to this subject have been published over the last three decades. There are several approximate methods to analyze vibration characteristics of these rotating beams. The finite element method (FEM) has been found

⇑ Corresponding author. E-mail address: [email protected] (Y.H. Li). http://dx.doi.org/10.1016/j.ymssp.2016.12.041 0888-3270/Ó 2017 Elsevier Ltd. All rights reserved.

24

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

in a large number of literature. Sai Ram and Sinha [5] studied the effects of temperature and moisture on the vibration characteristics of a laminated composite plate using FEM, in which the influence of transverse shear deformation was taken into account. According to the shear flexible theory, the free flexural vibration of a laminated composite rotating beam was investigated using FEM in Ref. [6]. Numerical results were presented, in which the modular ratio, rotation speed and slenderness ratio were taken into account. Based on the first order shear deformation theory, Parhi et al. [7] analyzed the free vibration and transient responses of a multiple delaminated composite shell in hygrothermal environment by a finite element formulation. Rao and Gupta [8] applied the finite element technique to obtain the natural frequencies and mode shapes in the bending-bending vibration of rotating Timoshenko beams. The effects of twist angle, shear deformation and rotation speed on the natural frequencies of the vibration were studied. A finite element analysis was presented for rotating cantilever beams in Ref. [9], and the linear partial differential governing equations of the rotating beam were derived using Hamilton’s principle. The effects of the rotation speed on the vibration of the beam were investigated. Naidu and Sinha [10] investigated the nonlinear free vibration of a laminated composite shell in hygrothermal environment by FEM. Maqueda et al. [11] obtained the eigenvalue solutions of a rotating beam by different nonlinear finite element formulations. The effect of centrifugal forces on these solutions was examined through numerical results. Liu and Jiang [12] studied the crack modeling of a rotating blade by an improved cracked finite element method. In addition, there exist other methods used by investigators for studying the vibrations of such beams. For example, the frequency and mode shape analysis of a rotating beam were presented by a semi-analytic method in papers [13,14]. Murthy [15] obtained the coupled flapwise bending and torsion vibration characteristics of a rotor blade adopting the integrating matrix method. Ӧzdemir and Kaya [16] acquired the natural frequencies of the flapwise bending vibration for a rotating cantilever Bernoulli-Euler beam using the differential transform method. The effects of some parameters such as angular velocity and hub radius on the vibration characteristics were discussed. Lang and Nemat-Nasser [17] determined the coupled vibration characteristics of a non-uniform rotating blade with the help of the transmission matrix method. Li et al. [18] solved the linear characteristic problem of a wind turbine blade applying the assumed modes method. And they discussed the effects of the rotation speed, pitch angle, and the wind velocity on flapwise dynamic response utilizing the multiplescales method. Surace et al. [19] analyzed the modal characteristics of a non-uniform rotating blade subjected to the coupled vibration by a numerical method based on Green’s function. Further, the effects due to hygrothermal environment on the vibration of the rotating beam were investigated by some methods in the following papers. Shen [20] investigated the effects of hygrothermal environment on postbuckling behaviors of a laminated plate by Reddy’s higher-order plate theory. Patel et al. [21] studied the static and dynamic characteristics of a thick laminated composite plate subjected to hygrothermal environment using a higher order theory. The effects of the temperature and moisture concentration on the hygrothermal behavior of a laminated composite plate were analyzed by a refined plate theory in Ref. [22]. Based on the D’Alembert’s principle, the governing equation of flapwise vibration for a rotating laminated composite beam under hygrothermal environment was obtained by Jiang in Ref. [23]. And the modal problem was solved employing the assumed-mode method. However, the influences of the design parameters such as the setting angle and the pitch angle on the vibration characteristics, which are essential in the engineering application of the wind turbine blade, were not considered. In this paper, the effects of the design parameters and hygrothermal environment will be all taken into account. The analysis of the vibration characteristics corresponding to the flapwise vibration of the rotating laminated composite beam is addressed. First, based on D’Alembert’s principle and the constitutive equation of one single layer material due to the hygrothermal effects, the governing equation of the rotating laminated composite beam in hygrothermal environment is established in Section 2, in which the setting angle and the pitch angle are considered. Then, a numerical approach based on Green’s functions is presented in Section 3 to obtain the vibration characteristics of the rotating beam. At last, in Section 4, the effects of the rotating speed, hub radius, fiber orientation, design parameters and hygrothermal environment on the natural frequencies are discussed.

2. Governing equation 2.1. Model of rotating composite beam Fig. 1 shows a schematic diagram of a rotating composite beam, attached to a rotating rigid hub radius R. In addition, the length, width and total thickness are represented by L, b and h respectively. bp denotes the pitch angle from the elastic axis to the vertical plane and h is the setting angle between rotating direction and the chord line. Several coordinate systems are used to describe the deformation of the beam. (xyz) is a body coordinate system, in which the x-axis corresponds to the undeformed elastic axis of the beam, the y-axis is along the lead-lag direction, and the z-axis lies along the flapwise direction. (XYZ) is a global coordinate system with the origin at the center of the volume of the hub. The X-axis is aligned with the horizontal direction. The Y-axis is parallel to the y-axis, and the Z-axis is perpendicular to the X-Y plane. Z-axis is the spindle of the rotating beam and K is the unit vector of it. The angular velocity can be expressed as X = XK, where X is the value of the rotating velocity. X will points to the direction of positive axis Z when the value of X is positive, otherwise direction of negative axis Z. Another new system (XrYrZr) can be obtained by translating the system (XYZ) from the hub center to the centroid of the fixed end of the beam, which will coincide with the system (xyz) if one rotate the Xr-Zr plane counterclockwise about the

25

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

Fig. 1. Schematic diagrams of a rotating composite cantilever beam. (a) Configuration of the rotating beam; (b) Pitch angle bp; (c) Setting angle h.

Yr-axis by the pitch angle bp. Hence the relationship between two unit vectors (i, j, k) and (I, J, K) of the two systems (xyz) and (XYZ) is given as

2 3 2 cos bp i 6 7 6 0 ¼ j 4 5 4  sin bp k

3 I 6 7 1 0 54 J 7 5 0 cos bp K 0

sin bp

32

ð1Þ

Fig. 2 shows the deformation of the section for the beam can be considered as the continuous changes of axial extension and flapwise bending. u and w stand for the axial and flapwise displacements of x- and z- directions for the rotating beam respectively. The axial extension yields a translation of the cross-section system (ngf), which does not alter the unit vector of the section system. Another new system (xsyszs) is obtained by translating the system (xyz) from the centroid of the fixed end of the beam to the shear center of the cross-section. Through the flapwise bending, the beam cross-section spins about the ys-axis by deflection angle hf (tanhf = w0 and (.)0 represents the derivative with respect to space coordinate x). The system (xsyszs) will turn into the system (ngf). Hence the relationship between two unit vectors (i, j, k) and (in, jg, kf) of the two systems (xyz) and (ngf) is given as

Fig. 2. Deformation of the beam section.

26

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

2

in

2

3

1

0

0

32

02

1  w2

76 6 7 6 4 jg 5 ¼ 4 0 cos h sin h 54 0 kf 0  sin h cos h w0

32 3 i 76 7 1 0 54 j 5 02 k 0 1  w2 0

w0

ð2Þ

The cross-section of the laminated composite beam is shown in Fig. 3. The number of layers is n and fiber orientation lie symmetrically with respect to the mid-plane of composite beam. Thickness of the k-th layer is defined as hk = zk  zk-1 and qk is the density of the material. 2.2. Basic equation of single layer material with hydrothermal effect The constitutive equation of single layer material in the hygrothermal environment can be given as [23]

ðrx Þk ¼ Q k11 ex  rT  rH

ð3Þ

where (rx)k, r and r represent total stress, thermal stress and humid stress of the k-th layer, respectively. ex denotes the total strain. The thermal stress generated by the temperature variation is expressed as T

H

rT ¼ Q k11 eTx þ Q k12 eTy þ Q k16 cTxy

ð4Þ

where the thermal strain components of Eq. (4) are defined as following forms

8 2 T 2 > > < ex ¼ ða1 cos H þ a2 sin HÞDT eTy ¼ ða1 sin2 H þ a2 cos2 HÞDT > > : cT ¼ 2ða  a Þ sin H cos HDT xy

1

ð5Þ

2

where a1 and a2 represent the thermal expansion coefficients along fiber orientation and direction normal to it, respectively. DT is the value of temperature variation. And H is the fiber orientation angle defined as the angle from axis x to the fiber orientation. The humid strain caused by the moisture concentration is expressed as

rH ¼ Q k11 eHx þ Q k12 eHy þ Q k16 cHxy

ð6Þ

where the humid strain components of Eq. (6) are defined as following forms

Fig. 3. Composition of a multi-layered composite beam. (a) Laminated composite beam geometry and layer numbering; (b) Fiber orientation angle of the kth layer.

27

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

8 2 H 2 > > < ex ¼ ðb1 cos H þ b2 sin HÞDC 2 eHy ¼ ðb1 sin H þ b2 cos2 HÞDC > > : cH ¼ 2ðb  b Þ sin H cos HDC 1

xy

ð7Þ

2

where b1 and b2 represent the moisture expansion coefficients along fiber orientation and normal direction, respectively. DC is the moisture concentration denoted as

DC ¼

DM DM f þ DM m ¼ M Mf þ Mm

ð8Þ

where M = Mf + Mm and DM = DMf + DMm are the mass of composite material in dry environment and the increment of mass in humid environment, respectively. Mf and Mm represent the mass of fiber and matrix. DMf and DMm represent the mass increment of fiber and matrix, respectively. In addition, Q k11 ; Q k12 and Q k16 in Eqs. (3), (4) and (6) are the elements of the reduced stiffness matrix Q of k-th layer as following [23]

8 2 4 k 4 2 > > < Q 11 ¼ ½Q 11 cos H þ 2ðQ 12 þ 2Q 66 Þ sin H cos H þ Q 22 sin Hk 2 4 Q k12 ¼ ½ðQ 11 þ Q 12  4Q 66 Þ sin H cos2 H þ Q 12 ðsin H þ cos4 HÞk > > : k 3 Q 16 ¼ ½ðQ 11  Q 12  2Q 66 Þ sin H cos3 H þ ðQ 12  Q 22 þ 2Q 66 Þ sin H cos Hk

ð9Þ

In Eq. (9),

Q 11 ¼ E1 =ð1  m12 m21 Þ; Q 22 ¼ E2 =ð1  m12 m21 Þ; Q 12 ¼ m21 E2 =ð1  m12 m21 Þ; Q 66 ¼ G12 where E1 and E2 are elastic modulus of material in principal orientation, G12 is shear module. v12 and ratios. Introducing Eqs. (4) and (6) into the Eq. (3), the generalized constitutive law can be obtained as

ð10Þ

v21 are the Poisson’s

ðrx Þk ¼ Q k11  ex  Ak DT  Bk DC

ð11Þ

where

Ak ¼ ½Q 11 ða1 cos2 H þ a2 sin HÞ þ Q 12 ða1 sin H þ a2 cos2 HÞ þ 2Q 16 ða1  a2 Þ sin H cos Hk 2

2

2

2

Bk ¼ ½Q 11 ðb1 cos2 H þ b2 sin HÞ þ Q 12 ðb1 sin H þ b2 cos2 HÞ þ 2Q 16 ðb1  b2 Þ sin H cos Hk 2.3. Axial force and moment of the composite beam The expressions of the axial force and moment for any cross-section of the beam can be obtained in this section. Considering the deformation of the composite beam in the x-z plane alone and the discontinuity of the stress components between each layers, the internal axial force along x-direction should be obtained by the integrations performed layer by layer as shown in the following form

Z Nx ¼

h=2 h=2

Z

b=2

b=2

ðrx Þk dydz ¼

Z

h=2

h=2

Z

b=2 b=2

  Q k11  ex  Ak DT  Bk DC dydz

ð12Þ

The axial normal strain-displacement relationship satisfies the following linear form [24]

ex ¼ u0  yw00 sin h  zw00 cos h

ð13Þ

where u and w denote the displacement components in x- and z-directions on the mid-plane. (.)0 and (.)00 mean the 1st and 2nd order partial derivatives with respect to x, respectively. Substituting Eq. (13) into Eq. (12), one has

Nx ¼ F T  NT  NH

ð14Þ

The part of the axial force Nx due to the elastic deformation is defined as

FT ¼

n Z X k¼1

zk

zk1

Z

b 2

2b

Q k11  ex dgdf ¼ EA  u0  EB  w00 cos h

where

EA ¼ b

n n X bX Q k11 ðzk  zk1 Þ EB ¼ Q k ðz2  z2k1 Þ 2 k¼1 11 k k¼1

Similarly, the forces NT and NH due to temperature variation and moisture concentration are formulated as

ð15Þ

28

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40 n Z X

NT ¼

k¼1

k¼1

b 2

2b

zk1

n Z X

NH ¼

Z

zk

Z

zk

zk1

b 2

2b

Ak DTdgdf ¼ bDTC NT

ð16Þ

Bk DCdgdf ¼ bDCC NC

ð17Þ

In the equations above, n n X X Ak ðzk  zk1 Þ; C NC ¼ Bk ðzk  zk1 Þ

C NT ¼

k¼1

k¼1

The bending moments of the composite beam of the cross-section in system (ngf) can be obtained as

Z Mg ¼

h 2

Z

2h

Z Mf ¼

h 2

h2

b 2

2b

Z

b 2

2b

 1  M ðrx Þk fdgdf ¼ EB  u0  EI1  w00 cos h  b DTC M T þ DCC C 2

ð18Þ

ðrx Þk gdgdf ¼ EI2 w00 sin h

ð19Þ

In the equations above,

CM T ¼ EI1 ¼

n n X X     Ak z2k  z2k1 ; C M Bk z2k  z2k1 C ¼ k¼1 n X  b Q k11 z3k 3 k¼1

k¼1



3

b  z3k1 ; EI2 ¼ 12

n X Q k11 ðzk  zk1 Þ k¼1

Translating the bending moments Mg and Mf to the coordinate system (xyz) through the relationship of Eq. (2) and neglecting the high-order terms, it can be obtained as

2 32 3 2 3 1 0 w0 i i   M ¼ ½ Mx M y Mz 4 j 5 ¼ 0 Mg Mf 4 w0 sin h cos h sin h 54 j 5 w0 cos h  sin h cos h k k  2  3T M M 1 b DTC T þ DCC C sin hw0  EB  sin hu0 w0  ðEI2  EI1 Þ sin h cos hw0 w00 2 3 2 6 7 i   6 7 2 M cos h 7 4 j 5 ¼ 6 EB  u0 cos h  ðEI1 cos2 h þ EI2 sin hÞw00  2b DTC M T þ DCC C 4 5   k M EB  u0 sin h  12 b DTC M sin h þ ðEI2  EI1 Þ sin h cos hw00 T þ DCC C

ð20Þ

For an inextensible beam, the geometric relation u0  wx2/2 holds. From Eqs. (14)–(17) and (20), the following results can be derived

Nx ¼

  1 EA  w02  EB  w00 cos h  b DTC NT þ DCC NC 2

ð21Þ

My ¼

 1 b 2 M DTC M EB  w02 cos h  ðEI1 cos2 h þ EI2 sin hÞw00  cos h T þ DCC C 2 2

ð22Þ

The centrifugal inertia force due to rotational motion and its two components in x- and z- directions respectively defined as

Z P¼ x

L

 1 mbðR þ x cos bp ÞX2 dx ¼ mbX2 RðL  xÞ þ ðL2  x2 Þ cos bp 2

ð23:aÞ

 1 PðxÞ ¼ mbX2 RðL  xÞ cos bp þ ðL2  x2 Þ cos2 bp 2

ð23:bÞ

 1 PðzÞ ¼ mbX2 RðL  xÞ sin bp þ ðL2  x2 Þ sin 2bp 4

ð23:cÞ

P P P P
and A are defined as q
¼ nk¼1 qk hk = nk¼1 hk and A ¼ b nk¼1 hk , denoting the average density of the where m ¼ b nk¼1 qk hk . q beam and the area of the cross-section, respectively. So the total axial force in the cross-section of rotating beam is presented as

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

  1 N ¼ F T  b DTC NT þ DCC NC þ mX2 ðL2  x2 Þ cos2 bp þ mRX2 ðL  xÞ cos bp 2

29

ð24Þ

where FT term in the right of the equation above denotes the axial force parts owing to the tension and bending. The 2nd and 3rd terms represent the axial forces due to the temperature variation and moisture concentration, respectively. The 4th and 5th terms stand for the axial force resulting from the rotation motion and the combined effect of the rotation motion and the hub radius, respectively. 2.4. Governing equation The free vibration equation of a rotating laminated composite beam subjected to the axial force is given as 0

€ ¼0 M00y þ ðNw0 Þ  mw

ð25Þ

Substituting Eqs. (22) and (24) into Eq. (25), the governing equation for the composite beam can be obtained as

h  i 3 2 ðEI1 cos2 h þ EI2 sin hÞwiv  EA  w02 w00 cos h þ b DTC NT þ DCC NC w00 2  1 2 2 2 2 € ¼0  mX ðL  x Þ cos bp þ mRX2 ðL  xÞ cos bp w00 þ mX2 ðR þ x cos bp Þ cos bp w0 þ mw 2

ð26Þ

where EA is the axial extension stiffness. EI1 and EI2 are the flapwise bending rigidity and the lagwise bending rigidity, respectively. In Eq. (26), hygrothermal environment and installation mode are taken into account for a rotating laminated composite beam. If the setting angle and the pitch angle are all ignored, the model problem (26) will be the same with that given by Jiang et al. [23]. There is one nonlinear term 3=2EA  w02 w00 cos h in Eq. (26). 2.5. Boundary conditions The geometric boundary conditions of flapwise vibration (26) for the rotating composite beam at the fixed end can be written as

x¼0:

w ¼ 0; w0 ¼ 0

ð27Þ

According to Eq. (22) and the relationship between shear force and bending moment, the following force boundary conditions at the free end can be expressed as

  2 M x ¼ L : 1=2EB  w02 cos h  ðEI1 cos2 h þ EI2 sin hÞw00  b=2 DTC M cos h ¼ 0; T þ DCC C 2

EB  w0 w00 cos h  ðEI1 cos2 h þ EI2 sin hÞw000 ¼ 0

ð28Þ

2.6. Dimensionless For the subsequent results to be general, the following dimensionless variables and parameters are introduced


¼ w

w x R ;
x ¼ ; R ¼ ;
t ¼ t L L L

sffiffiffiffiffiffiffiffiffi EI mL4 ; X2 ¼ X2 4 EI mL

ð29Þ

where EI = EI1cos2h + EI2cos2h is the stiffness of the beam. As a result, the dimensionless form of Eq. (26) can be written as follows


iv  a1 w
02 w
00 þ a2 w
00  w

 1 2 €
¼ 0
00 þ X2 ðR þ
x cos bp Þ cos bp w
0 þ w X ð1 
x2 Þ cos2 bp þ RX2 ð1 
xÞ cos bp w 2

ð30Þ

where a1 and a2 stand for the nonlinear term and the hygrothermal term respectively, whose expressions are given as following

  b DTC NT þ DCC NC L2 3EA cos h ; a2 ¼ a1 ¼ 2EI EI

ð31Þ

The dimensionless form of the boundary conditions are given by

x¼0: x¼1:


¼ 0 and w
0 ¼ 0 w h

 i M
02 cos h  b=2 DTC M
00 ¼ 0; 1=2EB  w cos h  L  EIw T þ DCC C
0w
00 cos h  L  EIw
000 ¼ 0 EB  w

ð32Þ ð33Þ

30

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40



For convenience, the dimensionless symbols of w; x; R;
t and X will be replaced by the common symbols of w, x, R, t and X in the following sections. 3. Modal problem In order to analyze the flapwise vibration characteristics of the rotating composite beam, it can be obtained only concerning the linear parts and neglecting the nonlinear term of Eq. (30)



€ wiv þ w

1 2 X ð1  x2 Þ cos2 bp þ RX2 ð1  xÞ cos bp w00 þ a2 w00 þ X2 ðR þ x cos bp Þ cos bp w0 ¼ 0 2

ð34Þ

By assuming that w(x, t) = c(x)eixt, where c(x) represents the modal function of the beam, x denotes the corresponding natural frequency and i is the imaginary unit. One can derive the following modal problem from Eq. (34)

" # 2 2 2 d d cðxÞ dcðxÞ d cðxÞ 2 2 ¼  x c ðxÞ  X ðR þ x cos b Þ cos b a 2 p p dx dx2 dx2 dx2  2 1 d cðxÞ þ X2 ð1  x2 Þcos2 bp þ RX2 ð1  xÞ cos bp 2 dx2

ð35Þ

The solution of Eq. (35) can be written as the following integral form [19]

cðxÞ ¼

Z

1

Gðx; nÞpðnÞdn

ð36Þ

0

where

pðxÞ ¼ x2 cðxÞ  X2 ðR þ x cos bp Þ cos bp

 2 dcðxÞ 1 d cðxÞ  a2  X2 ð1  x2 Þ cos2 bp  RX2 ð1  xÞ cos bp dx 2 dx2

ð37Þ

and G(x, n) is the Green’s function meaning the deflection w(x) at x due to a unit force applied at an arbitrary n. The expression of the Green’s function is given as following

Z

minðx;nÞ

Gðx; nÞ ¼

ðx  n1 Þðn  n1 Þdn1

ð38Þ

0

Substituting Eqs. (37) and (38) into Eq. (36), the integro-differential form can be obtained as following

cðxÞ ¼ x2

Z

1

0

Gcdn  X2

Z

1

GðR þ n cos bp Þ cos bp 0

Z

Z

1

Ga2

0

d c 1 dn þ X2 2 dx2 2

Z

Gð1  n2 Þ cos2 bp 0

d c dn dx2 2

1

d c dn dx2 2

1

þ RX2

dc dn  dx

Gð1  nÞ cos bp 0

ð39Þ

By discrete variable x along the interval [0, L], the integral and differential terms in Eq. (39) can be expressed approximately as following

R1 0

 df ðxÞ dx 

f ðnÞdn ¼

S X f iW i; i¼1

x¼xi

¼

S X

ð1Þ

dij f j ;

j¼1



d f ðxÞ  dx2 x¼x 2

¼ i

S X ð2Þ dij f j : j¼1

h i ð1Þ where fi = f(xi) and Wi is the weighting factor depending on the numerical integration method. It noted that D1 ¼ dij

SS

D2 ¼

ð2Þ ½dij SS

and

which represent the differential operators respectively.

Hence, Eq. (39) becomes



1 2

c ¼ x2 GWc  GWQD2 c  X2 GW M1 D1  M2 D2 c þ RX2 GWM3 D2 c where

c ¼ ½cðxi ÞS1 ; W ¼ diagðW i ÞSS ; G ¼ ½Gðxi ; xj ÞSS Q ¼ a2  ISS ; M1 ¼ cos bp  diagðR þ xi cos bp ÞSS M2 ¼ cos2 bp  diagð1  x2i ÞSS ; M3 ¼ cos bp  diagð1  xi ÞSS

ð40Þ

31

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

In above, the matrix I is unit vector. The matrix Q stands for the term resulting from temperature variation and moisture concentration. The matrix M1 represents the term due to the combined effect of pitch angle and hub radius. The matrix M2 and M3 are the terms caused by pitch angle. Eq. (40) can be expressed as the following form, which stands for a standard eigenvalue problem

½A1 ðI þ B1 þ X2 B2  RX2 B3 Þ  x2 Ic ¼ 0

ð41Þ

where

A ¼ GW;

B1 ¼ GWQD2 ;



1 B2 ¼ GW M1 D1  M2 D2 ; 2

B3 ¼ GWM3 D2

The natural frequencies of the beam can be obtained through solving the Eq. (41). 4. Numerical results and discussions In this section, the effects of several parameters on the flapwise free vibration characteristics of a rotating composite beam are explored. With the purpose of enabling a general application of results, they are presented in non-dimensional form by adopting the following non-dimensional natural frequency and rotating speed parameters:

xi ¼

xi X ; X ¼ ;x ¼ x0 x0 0

sffiffiffiffiffiffiffiffiffi EI

ð42Þ

mL4

To obtain the numerical results the composite beam is considered for this investigation, which is composed of 4 symmetric fiber orientation layers (H/H)s. The thickness of every layer is taken into being identical. The composite beam is made of a graphite-epoxy and the variation of the elastic modulus of graphite-epoxy in different hygrothermal environment has been provided in Table 1 and the density of the beam qk = 1600 kg/m3 in analysis. If there are no special instructions, the fiber orientation angle is taken into account as H = 45° in the following subsections. In addition, the length-to-radius ratio of the rotating beam is L/R = 25, and the thickness-to-radius ratio h/R = 0.25. The weighting matrix W and the two differential matrices D1 and D2 are determined by the Simpson’s formula and the central-difference method respectively. 4.1. Convergence and validity verification of the method In this subsection, the convergence of the first four natural frequencies derived by the numerical method based on Green’s function is verified in Fig. 4, where X⁄ = 5, DT = 25 °C, DC = 1.0%, h = 0o and bp = 0°. S in the figure stands for the number of collocation points. For the numerical method based on Green’s function, the accuracy of numerical results depends on the number S of collocation points. One can find that the first four natural frequencies converge rapidly. It can be confirmed that S = 200 is suitable to obtain reasonable values for the first four natural frequencies. Therefore the numerical results in the following subsections are obtained in the condition that S equals to 200. In Table 2 the results of the first four natural frequencies derived by the numerical method based on Green’s function are compared with results by the finite element method (FEM). In the modeling process of FEM, the beam elements with two nodes are applied through the Nastran software and the total number is 800. It can be found that results obtained from these two methods are extremely close, which means that the model in this study is reasonable and the numerical method based on Green’s function is valid. In order to further verify this numerical method based on Green’s function. By degenerating the governing equation for the rotating beam, the first four non-dimensional natural frequencies are compared with the results obtained from Ref. [25]

Table 1 Elastic modulus of graphite/epoxy [21]. Modulus (GPa)

E1 E2 G12 Modulus (GPa)

E1 E2 G12

(a) Moisture concentration, DC (%) 0.00

0.25

0.50

0.75

1.00

1.25

1.50

130 9.5 6.0

130 9.25 6.0

130 9.0 6.0

130 8.75 6.0

130 8.5 6.0

130 8.5 6.0

130 8.5 6.0

(b) Temperature, T (°C) 25

50

75

100

125

150

130 9.5 6.0

130 9.0 6.0

130 8.5 5.5

130 8.0 5.0

130 7.5 4.75

130 7.5 4.5

(a) At different moisture concentration: G13 = G12, G23 = 0.5 G12, m12 = 0.3, b1 = 0, b2 = 0.44  106. (b) At different temperatures: G13 = G12, G23 = 0.5 G12, m12 = 0.3, a1 = 0.3  106/°C, a2 = 28.1  106/°C.

32

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

5.895

24.69

1st Natural Frequency

24.68

ω*2

ω*1

5.8945 5.894 5.8935 5.893

2nd Natural Frequency

24.67 24.66

0

200

400

600

S

800

24.65

1000

0

200

400

600

64.5

4th Natural Frequency

64.45

123.8

ω4*

ω*3

1000

123.9 3rd Natural Frequency

64.4 64.35 64.3

800

S

123.7 123.6

0

200

400

600

800

1000

123.5

0

200

400

S

600

800

1000

S

Fig. 4. Convergence of the first four natural frequencies. (X = 5, DT = 25 °C, DC = 1.0%, h = 0°, bp = 0°). *

Table 2 Comparison of results from two methods (X* = 5, DT = 25 °C, DC = 1.0%, h = 0°, bp = 0°, S = 200).

x* i

FEM

Present

Error

x* 1 x* 2 x* 3 x* 4

5.894 24.677 64.380 123.671

5.893 24.680 64.394 123.709

0.02% 0.02% 0.02% 0.03%

in Table 3 when X⁄ = 0. As can be seen, the results derived by the numerical method based on Green’s function are very close to those of the reference. 4.2. Effects of rotating speeds on natural frequencies Table 4 gives the first four natural frequencies of the rotating beam at four different rotating speeds, where DT = 25 °C, DC = 1.0%, h = 0° and bp = 0°. As shown in this table, the first four natural frequencies all increase as the rotating speed increases. It also can be seen that the influences of the rotating speed on low-order natural frequencies are remarkable while these influences on high-order natural frequencies are small. This trend is also obtained in Ref. [26]. 4.3. Effects of setting angle on natural frequencies In Fig. 5, for selected values of the setting angle it is displayed in succession the variation of the first four natural frequencies as the rotating speed, where bp = 0°, DT = 25 °C and DC = 1.0%. It can be clearly found that the trend of variation of the first four natural frequencies as the rotating speed is comparatively evident, which is consistent with the trend supplied in

Table 3 Comparison of values for the first four natural frequencies (X* = 0).

x* i

Present

Ref. [25]

x* 1 x* 2 x* 3 x* 4

3.515 22.033 61.693 120.890

3.516 22.034 61.697 –

33

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40 Table 4 First four frequencies corresponding to different rotating speeds (DT = 25 °C, DC = 1.0%, h = 0°, bp = 0°).

X*

x* 1

x* 2

x* 3

x* 4

0 5 10 15

3.515 5.893 9.969 14.269

22.033 24.680 31.279 39.863

61.693 64.394 71.822 82.577

120.890 123.709 131.739 143.990

80

*

*

ω1

14

100

θ=0° θ=30 ° θ=60 ° θ=90 °

ω2

18

10

6

θ=0° θ=30 ° θ=60 ° θ=90 °

60

40

2 0

5

10

20 0

15

5

Ω* 300

250

10

15

10

15

Ω* 600

θ=0° θ=30 ° θ=60 ° θ=90 °

500

400 *

*

ω4

ω3

200

θ=0° θ=30 ° θ=60 ° θ=90 °

150

300

100

200

50 0

5

Ω*

10

15

100 0

5

Ω*

Fig. 5. Effects of setting angle on natural frequencies as different rotating speeds. (bp = 0°, DT = 25 °C, DC = 1.0%).

Table 3. Fig. 6 displays the variation of the first four natural frequencies as the setting angle, where the rotating speed is taken to be X⁄ = 5. As it clearly appears from this figure, the first four natural frequencies of the rotating beam all increase as the setting angle increases. In addition, the y-axis of coordinate system (xyz) will be parallel to the Y-axis of global coordinate system (XYZ) when the setting angle h = 0°. In this case, the stiffness of the rotating beam is EI1. The y-axis will rotate counterclockwise about the x-axis. When the setting angle h = 90°, the y-axis will be parallel to the Z-axis. In this case, the stiffness of the rotating beam is EI2. 4.4. Effects of pitch angle on natural frequencies In this subsection, the effects of the pitch angle on the first four natural frequencies as rotating speeds are presented in Fig. 7, where h = 30°, DT = 25 °C and DC = 1.0%. It shows that the pitch angle has significant effects on the first four natural frequencies. The reason for these impacts is that the variation of the pitch angle can lead to the changing the centrifugal force. The pitch angle affects the natural frequencies depending on the rotating speed. There is no effect of the pitch angle on the natural frequencies when X⁄ = 0. However, the influences of the pitch angle become large as the rotating speed increases. Fig. 8 displays the variation of the first four natural frequencies as the pitch angle, where the rotating speed is taken to be X⁄ = 5. It can be clearly found that in the figure, the first four natural frequencies of the rotating beam all decrease as the pitch angle increases, which is quite different from the effects of the setting angle.

34

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

12

70 60

10 *

ω2

ω*1

50 8

40 6

30 1st Natural Frequency

4

0

30

θ (°)

60

2nd Natural Frequency

20

90

0

30

60

90

θ (°) 400

150

300

ω*4

ω*3

200

200

100

4th Natural Frequency

3rd Natural Frequency

50

0

30

60

θ (°)

100

90

0

30

θ (°)

60

90

Fig. 6. Variation of natural frequencies as setting angle. (X* = 5, bp = 0°, DT = 25 °C, DC = 1.0%).

16

52

βp =0° βp =30 °

βp =30 °

βp =45 °

48

βp =60 °

βp =45 °

*

ω1

*

βp =60 °

ω2

12

βp =0°

44

8 40

4 0

122

5

Ω*

10

36 0

15

224

βp =0° βp =30 °

220

10

15

10

15

βp =0° βp =45 °

*

βp =60 °

ω4

ω3

*

βp =60 ° 114

110

106

Ω*

βp =30 °

βp =45 °

118

5

216

212

0

5

*

Ω

10

15

208 0

5

Ω*

Fig. 7. Effects of pitch angle on natural frequencies as different rotating speeds. (h = 30°, DT = 25 °C, DC = 1.0%).

35

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

40

7.5

39.5

ω*1

ω*2

8

39

7

2nd Natural Frequency

1st Natural Frequency

6.5

0

15

30

45

38.5

60

0

15

30

βp (°)

βp (°)

108.5

211

ω4*

211.5

ω*3

109

108

60

45

60

210.5 210

107.5

4th Natural Frequency

3rd Natural Frequency

107

45

0

15

30

45

βp (°)

209.5

60

0

15

30

βp (°)

Fig. 8. Variation of natural frequencies as pitch angle. (X* = 5, h = 30°, DT = 25 °C, DC = 1.0%).

16

52

ΔT=0°C ΔT=25 °C ΔT=50 °C

14

48 46 *

ω2

*

ω1

12 10 8

44 42

6.7 6.6

40

6.5

6

6.4

38

6.3 2.2

4

ΔT=0°C ΔT=25 °C ΔT=50 °C

50

0

124

5

Ω*

2.3

10

2.4

2.5

36

15

5

10

15

10

15

Ω* 225

ΔT=0°C ΔT=25 °C ΔT=50 °C

120

0

ΔT=0°C ΔT=25 °C ΔT=50 °C

220

ω4

*

ω3

*

116 215

112 210

108

104

0

5

*

Ω

10

15

205

0

5

Ω*

Fig. 9. Effects of temperature variation on natural frequencies. (h = 30°, bp = 10°, DC = 0.0%).

36

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

16

52 ΔC=0.0 % ΔC=1.0 %

14

ΔC=0.0 % ΔC=1.0 %

50 48 46 *

ω2

ω1

*

12 10

44 42

8

40 6 4

38 0

5

10

*

36 0

15

5

Ω 122

225

ΔC=0.0 % ΔC=1.0 %

*

114

110

106

15

ΔC=0.0 % ΔC=1.0 %

220

ω4

ω3

*

118

10

Ω*

215

210

0

5

10

15

205

0

5

*

Ω

Ω*

10

15

Fig. 10. Effects of moisture concentration on natural frequencies. (h = 30°, bp = 10°, DT = 0 °C).

7.73

40.1 2nd Natural Frequency

7.72

40

7.71

39.9

ω2*

ω*1

1st Natural Frequency

7.7

39.8

7.69 0

0.5

39.7

1

ΔC (%)

1

213

3rd Natural Frequency

4th Natural Frequency

109.2

212.5

109

ω*4

ω*3

0.5

ΔC (%)

109.4

108.8 108.6

0

212

211.5

0

0.5

ΔC (%)

1

211

0

0.5

1

ΔC (%)

Fig. 11. Variation of natural frequencies with moisture concentration. (X* = 5, h = 30°, bp = 10°, DT = 0 °C).

37

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

4.5. Effects of temperature variation on natural frequency In Fig. 9 the influences of the temperature variation on the first four natural frequencies of the rotating beam as the rotating speed are illustrated, where h = 30°, bp = 10° and DC = 0.0%. Thermal environment can determine the elasticity and thermoplastic properties of structures. The changing of elastic moduli due to the temperature variation and thermal expansion deformation has been considered in this subsection. As depicted in Fig. 9, the first four natural frequencies reduce with the increase of the temperature variation. In addition, effects of the temperature variation on first natural frequencies are very little while these effects on high-order ones become dramatic. The temperature variation has relatively large influence on the 400 1st Natural Frequency 2nd Natural Frequency 3rd Natural Frequency 4th Natural Frequency

300

14 12

ω*i

10 8

200

6 4

0

30

60

90

100

0

0

30

60

o

Θ( )

90

Fig. 12. Variation of natural frequencies with fiber orientation angle. (X* = 5, h = 30°, bp = 10°, DT = 25 °C, DC = 1.0%).

12

45

11

44 43

Ω* = 0 Ω* = 5 Ω* = 10

10

2

ω*

1

42

ω*

9 8

39

6

38 0

0.02

0.04

R

0.06

0.08

37

0.1

114

0.02

0.04

R

0.06

215

4

Ω* = 0 Ω* = 5 Ω* = 10

ω*

3

ω*

0

0.08

0.1

217

112

110

108

106

41 40

7

5

Ω* = 0 Ω* = 5 Ω* = 10

Ω* = 0 Ω* = 5 Ω* = 10

213

211

0

0.02

0.04

R

0.06

0.08

0.1

209 0

0.02

0.04

R

0.06

0.08

Fig. 13. Variation of natural frequencies with hub radius. (h = 30°, bp = 10°, DT = 25 °C, DC = 1.0%).

0.1

38

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

high-order natural frequencies when the rotating speed is small. However, this influence of the temperature variation will narrow as the increase of the rotating speed. 4.6. Effects of moisture concentration on natural frequencies For selected values of the moisture concentration the variations of the first four natural frequencies are given in Fig. 10, where h = 30°, bp = 10° and DT = 0 °C. The results reveal that the differences between the natural frequencies corresponding to DC = 0.0% and DC = 1.0%. It clearly shows that moisture concentration has fairly small effects on the first four natural frequencies of the rotating beam. Fig. 11 demonstrates the effects of moisture concentration from 0% to 1.0% on the first four natural frequencies of the rotating beam in the case that X⁄ = 5. Comparing the influences of the temperature variation the moisture concentration has weak effects on the first four natural frequencies of the beam, especially on the low-order ones. 4.7. Effects of fiber orientation angle on natural frequencies Fig. 12 shows the variations of the first four natural frequencies as the fiber orientation angle increasing from 0° to 90° with fixed conditions of X⁄ = 5, h = 30°, bp = 10°, DT = 25 °C and DC = 1.0%. In this figure, the horizontal ordinate stands for fiber orientation angle of laminate, as for example ‘‘60” is (60°/60°)s. It is observed that the first four natural frequencies decrease as the fiber orientation increases. And the first four natural frequencies all reduce rapidly in the range of 30– 60°. However, their reduction rates become smaller when the fiber orientation angle exceed higher than 60°. 4.8. Effects of hub radius on natural frequencies For selected values of the rotating speed, Fig. 13 displays the variation of the first four natural frequencies of the rotating beam as the hub radius, where h = 30°, bp = 10°, DT = 25 °C and DC = 1.0%. The results reveal a continuous linear increase of the first four natural frequencies that accompanies the increase of the hub radius. The reason is that the centrifugal inertia force as a result of the hub radius strengthens the bending stiffness of the rotating beam, which yields the increase of the natural frequencies. In addition, the values of the first four natural frequencies have nothing to do with the hub radius when X⁄ = 0. The explanation for this phenomenon in Fig. 13 is that the hub radius effects wear off when X⁄ = 0, which can be also explained by the expression of Eq. (24) in Section 2. 5. Conclusions In the present paper, the modal problem for a rotating laminated composite beam under hygrothermal environment has been handled, in which the setting angle and pitch angle are considered. Combined effects of the rotating speed, the setting angle, the pitch angle, the temperature variation, the moisture concentration, the fiber orientation and the hub radius on vibration characteristics of flapwise vibration for the rotating composite beam are investigated by the numerical method based on Green’s function. The following conclusions can be drawn: (1) Natural frequencies of the rotating beam ascend as the increase of the rotating speed. Effects of the rotating speed on low-order frequencies are remarkable than those on high-order frequencies. (2) Influences of the setting angle and the pitch angle on natural frequencies are mostly by means of the rotating speed. The natural frequencies all increase as the setting angle increases, while decrease as the pitch angle increases. (3) The temperature variation has more significant effects on the natural frequencies than the moisture concentration does, especially on the high-order ones. (4) Influences of the fiber orientation angle and hub radius on the natural frequencies are all great.

Acknowledgments This research work was supported by the National Natural Science Foundation of China (Grant Nos. 11372257, 11302183 and 11102030). Appendix Details on the Eq. (13) After deformation, the position vector of an arbitrary point of the beam can be expressed as:

2 3 2 3 i i 6 7 6 7 R ¼ ½ x þ u; v ; w 4 j 5 þ ½ 0; y; z T4 j 5 k k

ðA:1Þ

39

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

where T is the transformation matrix of the relationship between two unit vectors (i, j, k) and (in, jg, kf) of the two systems (xyz) and (ngf), whose expression is given as following

2

32 02 1 0 0 1  w2 6 76 T ¼ 4 0 cos h sin h 54 0 0  sin h cos h w0

0

w0

1

0

3 7 5

ðA:2Þ

02

0 1  w2

Substituting Eq. (A.2) into Eq. (A.1) and ignoring the non-linear term, it yields

R ¼ ðx þ u  y sin hw0  z cos hw0 Þi þ ðv þ y cos h  z sin hÞj þ ðw þ y sin h þ z cos hÞk

ðA:3Þ

Differentiating the Eq. (A.3), it can be obtained

dR ¼ ð1 þ u0  y sin hw00  z cos hw00 Þi þ v 0 j þ w0 k dx

ðA:4Þ

The deformation of the beam depends on the axial strain mainly, whose can be derived as

e¼

d
s=dx  ds=dx ds=dx

ðA:5Þ

where d
s/dx stands for the length of dR/dx.

ð1 þ xÞa ¼ 1 þ ax þ

aða  1Þ 2!

x2 þ

aða  1Þða  2Þ 3!

x3 þ   

ðA:6Þ

Using the approximate expression (A.6) and ignoring the higher order infinitesimal terms,

  d
s dR v 02 w02 ¼   ¼ 1 þ u0 þ þ  y sin hw00  z cos hw00 dx dx 2 2

ðA:7Þ

For the un-deformed ds/dx, it can be obtained when u = v = w = 0 and u = v = w = v = w = 0 0

ds ¼1 dx

0

0

00

00

ðA:8Þ

Hence,

e¼

d
s=dx  ds=dx v 02 w02 ¼ u0 þ þ  y sin hw00  z cos hw00 ds=dx 2 2

ðA:9Þ

Just considering the axial displacement u and flapwise displacement w and neglecting the non-linear terms, so it fields

e¼

d
s=dx  ds=dx ¼ u0  y sin hw00  z cos hw00 ds=dx

ðA:10Þ

References [1] E.L. Oliveira, N.M.M. Maia, A.G. Marto, R.G.A. Da Silva, F.J. Afonso, A. Suleman, Modal characterization of composite flat plate models using piezoelectrictransducer, Mech. Syst. Signal Process. 79 (2016) 16–29. [2] P. Kudela, W. Ostachowicz, A. Zak, Damage detection in composite plates with embedded PZT transducers, Mech. Syst. Signal Process. 22 (2008) 1327– 1335. [3] Y. Qin, X. Li, E.C. Yang, Y.H. Li, Flapwise free vibration characteristics of a rotating composite thinwalled beam under aerodynamic force and hygrothermal environment, Compos. Struct. 153 (2016) 490–503. [4] Ӧ. Ӧzdemir, M.O. Kaya, Energy derivation and extension-flapwise bending vibration analysis of a rotating piezolaminated composite Timoshenko beam, Mech. Adv. Mater. Struct. 21 (6) (2014) 477–489. [5] K.S. Sai Ram, P.K. Sinha, Hygrothermal effects on the bending characteristics of laminated composite plates, Comput. Struct. 40 (4) (1991) 1009–1015. [6] B.P. Patel, M. Ganapathi, M. Touratier, Free vibrations analysis of laminated composite rotating beam using C1 shear flexible element, Defence Sci. J. 49 (1) (1999) 3–8. [7] P.K. Parhi, S.K. Bhattacharyya, P.K. Sinha, Hygrothermal effects on the dynamic behavior of multiple delaminated composite plates and shells, J. Sound Vib. 248 (2) (2001) 195–214. [8] S.S. Rao, R.S. Gupta, Finite element vibration analysis of rotating Timoshenko beams, J. Sound Vib. 242 (1) (2001) 103–124. [9] J. Chung, H.H. Yoo, Dynamic analysis of a rotating cantilever beam by using the finite element method, J. Sound Vib. 249 (1) (2002) 147–164. [10] N.V.S. Naidu, P.K. Sinha, Nonlinear free vibration analysis of laminated composite shells in hygrothermal environments, Compos. Struct. 77 (4) (2007) 475–483. [11] L.G. Maqueda, O.A. Bauchau, A.A. Shabana, Effect of the centrifugal forces on the finite element eigenvalue solution of a rotating blade: a comparative study, Multibody Syst. Dyn. 19 (3) (2008) 281–302. [12] C. Liu, D.X. Jiang, Crack modeling of rotating blades with cracked hexahedral finite element method, Mech. Syst. Signal Process. 46 (2014) 406–423. [13] R.O. Stafford, V. Giurgiutiu, Semi-analytic methods for rotating Timoshenko beams, Int. J. Mech. Sci. 17 (11) (1975) 719–727. [14] V. Giurgiutiu, R.O. Stafford, Semi-analytic methods for frequencies and mode shapes of rotor blades, Vertica 1 (4) (1977) 291–306. [15] V.R. Murthy, Dynamic characteristics of rotor blades: integrating matrix method, AIAA J. 15 (4) (1977) 595–597. [16] Ӧ. Ӧzdemir, M.O. Kaya, Flapwise bending vibration analysis of a rotating tapered cantilever Bernoulli-Euler beam by differential transform method, J. Sound Vib. 289 (1) (2006) 413–420. [17] K.W. Lang, S. Nemat-Nasser, An approach for estimating vibration characteristics of rotor blades, AIAA J. 17 (9) (1979) 995–1002.

40

Y. Qin, Y.H. Li / Mechanical Systems and Signal Processing 91 (2017) 23–40

[18] L. Li, Y.H. Li, H.W. Lv, Q.K. Liu, Flapwise dynamic response of a wind turbine blade in super-harmonic resonance, J. Sound Vib. 331 (17) (2012) 4025– 4044. [19] G. Surace, V. Anghel, C. Mares, Coupled bending-bending-torsion vibration analysis of rotating pretwisted blades: an integral formulation and numerical examples, J. Sound Vib. 206 (4) (1997) 473–486. [20] H.S. Shen, Hygrothermal effects on the postbuckling of shear deformable laminated plates, Int. J. Mech. Sci. 43 (5) (2001) 1259–1281. [21] B.P. Patel, M. Ganapathi, D.P. Makhecha, Hygrothermal effects on the structural behavior of thick composite laminates using higher-order theory, Compos. Struct. 56 (1) (2002) 25–34. [22] A.M. Zenkour, M.N.M. Allam, A.F. Radwan, Effects of hygrothermal conditions on cross-ply laminated plates resting on elastic foundations, Arch. Civ. Mech. Eng. 14 (1) (2014) 144–159. [23] B.K. Jiang, J. Xu, Y.H. Li, Flapwise vibration analysis of a rotating composite beam under hygrothermal environment, Compos. Struct. 117 (2014) 201– 211. [24] B.L. Li, X.G. Song, D.X. He, Y.H. An, Structural Dynamics of Wind Turbines, The BUAA Press, Beijing, 1999 (In Chinese). [25] J.R. Banerjee, H. Su, Development of a dynamic stiffness matrix for free vibration analysis of spinning beams, Comput. Struct. 82 (2004) 2189–2197. [26] Y.H. Li, L. Li, Q.K. Liu, H.W. Lv, Dynamic characteristics of lag vibration of a wind turbine blade, Acta Mech. Solida Sin. 26 (6) (2013) 592–602.