# A Corotational Finite Element Method Combined with Floating Frame Method for Large Steady-State Deformation and Free Vibration Analysis of a Rotating-Inclined Beam

Authors: Ming Hsu Tsai; Wen Yi Lin; Yu Chun Zhou; Kuo Mo Hsiao
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

A corotational finite element method combined with floating frame method and a numerical procedure is proposed to investigate large steady-state deformation and infinitesimal-free vibrationaround the steady-state deformation of a rotating-inclined Euler beam at constant angular velocity. The element nodal forces are derived using the consistent second-order linearization of the nonlinear beam theory, the d'Alembert principle, and the virtual work principle in a current inertia element coordinates, which is coincident with a rotating element coordinate system constructed at the current configuration of the beam element. The governing equations for linear vibration are obtained by the first-order Taylor series expansion of the equation of motion at the position of steady-state deformation. Numerical examples are studied to demonstrate the accuracy and efficiency of the proposed method and to investigate the steady-state deformation and natural frequency of the rotating beam with different inclined angle, angular velocities, radius of the hub, and slenderness ratios.


---

## Body

## 1. Introduction

Rotating beams are often used as a simple model for propellers, turbine blades, and satellite booms. Rotating beam differs from a nonrotating beam in having additional centrifugal force and Coriolis effects on its dynamics. The vibration analysis of rotating beams has been extensively studied [1–25].  However, the vibration analysis of rotating beam with inclination angle, which is considered in the recent computer cooling fan design on the natural frequencies of rotating beams [21], is rather rare in the literature [10, 19, 21, 22]. In [21, 22], the effect of the steady-state axial deformation and the inclination angle on the natural frequencies of the rotating beam was investigated. However, the lateral steady-state deformation and its effects on the natural frequencies of the rotating beam were not considered in [21, 22]. To the authors’ knowledge, the lateral steady-state deformation and its effects on the lagwise bending and axial vibration of rotating inclined beams are not reported in the literature.It is well known that the spinning elastic bodies sustain a steady-state deformation (time-independent deformation) induced by constant rotation [26].  For rotating beams with an inclination angle as shown in Figure 1, the steady-state deformations include axial deformation and lateral deformation.The linear solution of the steady-state deformation of rotating-inclined beam induced by constant rotation can be easily obtained using mechanics of materials.  However, the centrifugal stiffening effect on the steady lateral deformation is significant for slender rotating-inclined beam, and the centrifugal force is configuration dependent load; thus the linear solution of the steady-state deformation of rotating inclined beam may be not accurate enough. The lagwise bending and axial vibration of rotating inclined beams are coupled due to the Coriolis effects [15, 24] and the lateral steady-state deformation. The accuracy of the frequencies obtained from linearizing about the steady-state deformation is dependent on the accuracy of the steady-state deformation and the accuracy of the linearized perturbation [6, 12]. Thus, the geometrical nonlinearities that arise due to steady-state deformation should be considered. In [6], the rotating beam with pretwist, precone, and setting angle is studied. The undeformed state of the rotating beam is chosen to be the reference state to define the deformation parameters of the rotating beam. The geometric nonlinearities up to the second degree are considered. The Galerkin method, with vibration modes of nonrotating beam, is employed for the solution of both steady-state nonlinear equations and linear perturbation equations. In [8], it is reported that for a cantilever beam with a tip mass, even up to the third degree geometric nonlinearities are considered, in some cases, very inaccurate eigenvalues for the perturbed linearized equation of motion are obtained. The formulation used in [6, 8, 12] may be regarded as a total Lagrangian (TL) formulation combined with the floating frame method. In order to capture correctly all inertia effects and coupling among bending, twisting, and stretching deformations of the rotating beam, the governing equations of the rotating beam might be derived by the fully geometrically nonlinear beam theory [12, 27, 28]. The exact expressions for the inertia, deformation forces, and the governing equations of the rotating beam, which are required in a TL formulation for large displacement/small strain problems, are highly nonlinear functions of deformation parameters. However, the dominant factors in the geometrical nonlinearities of beam structures are attributable to finite rotations, with the strains remaining small. For a beam structures discretized by finite elements, this implies that the motion of the individual elements to a large extent will consist of rigid body motion. If the rigid body motion part is eliminated from the total displacements and the element size is properly chosen, the deformational part of the motion is always small relative to the local element axes; thus in conjunction with the corotational formulation, the higher-order terms of nodal deformation parameters in the element deformation and inertia nodal forces may be neglected by consistent linearization [28, 29]. In [29], Hsiao et al. presented a corotational finite element formulation and numerical procedure for the dynamic analysis of planar beam structures. Both the element deformation and inertia forces are systematically derived by consistent linearization of the fully geometrically nonlinear beam theory using the d’Alembert principle and the virtual work principle.  This formulation and numerical procedure were proven to be very effective by numerical examples studied in [29]. However, because the nodal displacements and rotations, velocities, accelerations, and the equations of motion of the system are defined in terms of a fixed global coordinate system, the formulation proposed in [29] cannot be used for steady-state deformation and free vibration analysis of a rotating-inclined beam. The absolute nodal coordinate formulation [30, 31] is used to large rotation and large deformation problems. Numerical results show that the absolute nodal coordinate formulation can be effectively used in the large deformation problems. However, the mass matrix of the finite elements in [30, 31] is a constant matrix, and therefore, the centrifugal and Coriolis forces are equal to zero. Thus, the absolute nodal coordinate formulation cannot be used for steady-state deformation and free vibration analysis of a rotating inclined beam.A rotating inclined beam, (a) top view, (b) side view.
(a)(b)The objective of this study is to present a corotational finite element method combined with floating frame method and a numerical procedure for large steady-state deformation and free vibration analysis of a rotating-inclined beam at constant angular velocity. The nodal coordinates, displacements and rotations, absolute velocities, absolute accelerations, and the equations of motion of the system are defined in terms of an inertia global coordinate system which is coincident with a rotating global coordinate system rigidly tied to the rotating hub, while the total deformations in the beam element are measured in an inertia element coordinate system which is coincident with a rotating element coordinate system constructed at the current configuration of the beam element. The rotating element coordinates rotate about the hub axis at the angular speed of the hub. The inertia nodal forces and deformation nodal forces of the beam element are systematically derived by the virtual work principle, the d’Alembert principle, and consistent second-order linearization of the fully geometrically nonlinear beam theory [27–29] in the element coordinates. Due to the consideration of the exact kinematics of Euler beam, some coupling terms of axial and flexural deformations are retained in the element internal nodal forces. The element equations are constructed first in the inertia element coordinate system and then transformed to the inertia global coordinate system using standard procedure.The infinitesimal-free vibrations of rotating beam are measured from the position of the corresponding steady-state deformation. The governing equations for linear vibration of rotating beam are obtained by the first-order Taylor series expansion of the equation of motion at the position of steady-state deformation.Dimensionless numerical examples are studied to demonstrate the accuracy and efficiency of the proposed method and to investigate the effect of inclination angle and slenderness ratio on the steady-state deformation and the natural frequency for rotating inclined Euler beams at different angular speeds.


## 2. Formulation

### 2.1. Description of Problem

Consider an inclined uniform Euler beam of lengthLT rigidly mounted with an inclination angle α on the periphery of rigid hub with radius R rotating about its axis fixed in space at a constant angular speed Ω as shown in Figure 1. The axis of the rotating hub is perpendicular to one of the principal directions of the cross section of the beam. The deformation displacements of the beam are defined in an inertia rectangular Cartesian coordinate system which is coincident with a rotating rectangular Cartesian coordinate system rigidly tied to the hub.Here only axial and lagwise bending vibrations are considered. It is well known that the beam sustains a steady-state deformations (time-independent deformation displacements) induced by constant rotation [26]. In this study, large displacement and rotation with small strain are considered in the steady-state deformation. The vibration (time-dependent deformation displacements) of the beam is measured from the position of the steady-state deformation, and only infinitesimal-free vibration is considered. Note that the axial and lagwise vibrations, which are coupled due to the Coriolis effects and the lateral steady-state deformation, cannot be analyzed independently. Here the engineering strain and stress are used for the measure of the strain and stress.


### 2.2. Basic Assumptions

The following assumptions are made in derivation of the beam element behavior.(1)
The beam is prismatic and slender, and the Euler-Bernoulli hypothesis is valid.(2)
The unit extension of the centroid axis of the beam element is uniform.(3)
The deformation displacements and rotations of the beam element are small.(4)
The strains of the beam element are small.In conjunction with the corotational formulation and rotating frame method, the third assumption can always be satisfied if the element size is properly chosen. Thus, only the terms up to the second order of deformation parameters and their spatial derivatives are retained in element position vector, strain, and deformation nodal forces by consistent second-order linearization in this study.


### 2.3. Coordinate Systems

In order to describe the system, we define three sets of right-handed rectangular Cartesian coordinate systems.(1)
A rotating global set of coordinates,Xi(i=1,2,3) (see Figures 1 and 2); the coordinates rotate about the hub axis at a constant angular speed Ω as shown in Figure 1. The origin of this coordinate system is chosen to be the intersection of the centroid axes of the hub and the undeformed beam. The X1 axis is chosen to coincide with the centroid axis of the undeformed beam, and the X2 and X3 axes are chosen to be the principal directions of the cross section of the beam at the undeformed state. The direction of the axis of the rotating hub is parallel to the X3 axis. The nodal coordinates, nodal deformation displacements, absolute nodal velocity, absolute nodal acceleration, and equations of motion of the system are defined in terms of an inertia global coordinate system which is coincident with the rotating global coordinate system.(2)
Element coordinates;xi(i=1,2,3) (see Figure 2), a set of element coordinates is associated with each element, which is constructed at the current configuration of the beam element. The coordinates rotate about the hub axis at a constant angular speed Ω. The origin of this coordinate system is located at the element node 1, the centroid of the end section. The  x1 axis is chosen to pass through two end nodes of the element; the directions of the x2 and x3 axes are chosen to coincide with the principal direction of the cross section in the undeformed state. Because only the displacements in X1X2 plane are considered, the directions of x3 axis and X3 axis are coincident. The position vector, deformations, absolute velocity, absolute acceleration, internal nodal forces, stiffness matrices, and inertia matrices of the elements are defined in terms of an inertia element coordinate system which is coincident with the rotating element coordinate system.Figure 2
Coordinate systems.In this study, the direction of the axis of the rotating hub is parallel to theX3 axis and only the displacements in X1X2 plane are considered. Thus, the angular velocity of the hub referred to the global coordinates may be given by(2.1)ΩG={00Ω},
where the symbol {} denotes a column matrix, which is used through the paper.


### 2.4. Kinematics of Beam Element

LetQ (Figure 3) be an arbitrary point in the beam element and P the point corresponding to Q on the centroid axis. The position vector of point Q in the undeformed configurations referred to the current element coordinate system may be expressed as(2.2)r0={x,y,z}.Figure 3
Kinematics of Euler beam.Using the approximationcosθ≈1-(1/2)θ2, sinθ≈θ, and (1+εc)≈1, retaining all terms up to the second order, the position vector of point Q in the deformed configurations referred to the current element coordinate system may be expressed as(2.3)r={r1,r2,r3}={xp-yθ,y(1-12θ2)+v,z},(2.4)θ≈sin⁡θ=∂v(x,t)∂s=∂v(x,t)∂x∂x∂s=v′1+εc≈v′,(2.5)εc=∂s∂x-1,
where xp(x,t) and v(x,t) are the x1 and x2 coordinates of point P, respectively, in the deformed configuration, t is time, θ=θ(x,t) is the angle counterclockwise measured from x1 axis to the tangent of the centroid axis of the deformed beam, εc is the unit extension of the centroid axis, and s is the arc length of the deformed centroid axis measured from node 1 to point P. In this paper, ()′ denotes (),x=∂()/∂x.Here, the lateral deflection of the centroid axis,v(x,t), is assumed to be the Hermitian polynomials of x and may be expressed by(2.6)v(x,t)={N1,N2,N3,N4}t{v1,v1′,v2,v2′}=Nbtub,
where vj=vj(t) and vj′=vj′(t)(j=1,2) are nodal values of v and v,x, respectively, at nodes j. Note that, due to the definition of the element coordinates, the values of vj(j=1,2) are zero. However, their variations and time derivatives are not zero. Ni(i=1–4) are shape functions and are given by(2.7)N1=14(1-ξ)2(2+ξ),N2=L8(1-ξ2)(1-ξ),N3=14(1+ξ)2(2-ξ),N4=L8(-1+ξ2)(1+ξ),(2.8)ξ=-1+2xL,
where L is the length of the undeformed beam element.Making use of assumptionsv,x≪1 and εc≪1, the relationship between xp(x,t), v(x,t), and x in (2.3) may be approximated by(2.9)xp(x,t)=u1+∫0x(1+εc-12v,x2)dx,
where u1 is the displacement of node 1 in the x1 direction. Note that due to the definition of the element coordinate system, the value of u1 is equal to zero. However, the variation and time derivatives of u1 are not zero.The axial displacements of the centroid axis may be determined from the lateral deflections and the unit extension of the centroid axis using (2.9).From (2.9), one may obtain
(2.10)l=L+u2-u1=xc(L,t)-xc(0,t)=∫0L(1+εc-12v,x2)dx
in which ℓ is the current chord length of the centroid axis of the beam element and u2 is the displacement of node 2 in the x1 direction. Using the assumption of uniform extension of the centroid axis and (2.10), εc in (2.10) maybe expressed by(2.11)εc=1L(Gatua+12Gbtub),Ga={-1,1},ua={u1,u2},Gb=∫0LNb′v,xdx.Substituting (2.11) into (2.9), one may obtain(2.12)xp(x,t)=Natua+x+x2LGbtub-12∫0xv,x2dx,Na={1-ξ2,1+ξ2}.From (2.3) and the definition of engineering strain [32, 33], making use of the assumption of small strain, and retaining the terms up to the second order of deformation parameters, the engineering strain in the Euler beam may be approximated by(2.13)ε11=εc-yv′′.The absolute velocity and acceleration vectors of pointQ in the beam element may be expressed as(2.14)v={v1,v2,v3}=vo+Ω×r+ṙ,(2.15)a={a1,a2,a3}=ao+Ω̇×r+Ω×(Ω×r)+2Ω×ṙ+r̈,(2.16)vo=Ω×rAo,(2.17)ao={ao1,ao2,ao3}=Ω×(Ω×rAo),(2.18)Ω=AGEtΩG,(2.19)rAo=AGEtrAoG,(2.20)rAoG=rAO+rOoG={Rcos⁡α+Xo,-Rsinα+Yo,0},
where r is the position vector of point Q given in (2.3) referred to the current moving element coordinate system, the symbol ()̇ denotes time derivative, Ω is the vector of angular velocity referred to the current inertia element coordinates, ΩG is the angular velocity of the hub referred to the global coordinates given in (2.1), AGE is the transformation matrix between the current global coordinates and the current element coordinates, vo and ao are the absolute velocity and absolute acceleration of point o, the origin of the current element coordinates, Xo and Yo are coordinates of point o referred to the current global coordinates, R is the radius of the hub, and α is inclination angle of the rotating beam.  Ω×(Ω×r)and 2Ω×ṙ are centripetal acceleration and Coriolis acceleration, respectively. ṙ and r̈ are the velocity and acceleration of point Q relative to the current moving element coordinates. From (2.3), (2.11) and (2.12), ṙ and r̈ may be expressed as
(2.21)ṙ={ṙ1,ṙ2,ṙ3}={ẋp-yv̇,x,v̇-yv̇,xv,x,0},r̈={r̈1,r̈2,r̈3}={ẍp-yv̈,x,v̈-yv̇,x2-yv̈,xv,x,0},ẋp=Natu̇a+xLGbtu̇b-∫0xv,xv̇,xdx,ε̇c=1L(Gatu̇a+Gbtu̇b),ẍp=Natüa+xL(Gbtüb+Ġbtu̇b)-∫0x(v,xv̈,x+v̇,xv̇,x)dx,ε̈c=1L(Gatüa+Ġbtu̇b+Gbtüb).Note that the current element coordinates constructed at the current configuration of the beam element rotate about the hub axis at the angular velocity of the hub. Thus, the centripetal acceleration and Coriolis acceleration corresponding to the inertia forces of the rotating beam are unique. For nonrotating beam,Ω=0 and ṙ and r̈ are the absolute velocity and acceleration referred to the current element coordinate.


### 2.5. Element Nodal Force Vector

Letuj,δvj, and δvj′(j=1,2) denote the virtual displacements in the x1 and x2 directions of the current inertia element coordinates, and virtual rotations applied at the element nodes j. The element nodal force corresponding to virtual nodal displacements δuj, δvj, and δvj′(j=1,2) are fij, the forces in the xi(i=1,2) directions, and mj moments about the x3 axis, at element local nodes j.The element nodal force vector is obtained from the d’Alembert principle and the virtual work principle in the current inertia element coordinates. The virtual work principle requires that(2.22)δuatfa+δubtfb=∫V(δε11σ11+ρδrta)dV,(2.23)δua={δu1,δu2},(2.24)δub={δv1,δv1′,δv2,δv2′},(2.25)fa=faD+faI={f11,f12},(2.26)fb=fbD+fbI={f21,m1,f22,m2},(2.27)faD={f11D,f12D},(2.28)fbD={f21D,m1D,f22D,m2D},(2.29)faI,={f11I,f22I},(2.30)fbI,={f21I,m1I,f22I,m2I},
where fi(i=a,b) are the generalized force vectors corresponding to δua and δub, respectively, fiD and fiI(i=a,b) are element deformation nodal force vector and inertia nodal force vector corresponding to fi, respectively, V is the volume of the undeformed beam element, and δε11 is the variation of ε11 in (2.13) corresponding to δua and δub. σ11 is the engineering stress. For linear elastic material, σ11=Eε11, where E is Young’s modulus. ρ is the density, δr is the variation of r in (2.3) (referred to the current inertia element coordinate system) corresponding to δua and δub, and a is the absolute acceleration in (2.15).If the element size is chosen to be sufficiently small, the values of the deformation parameters of the deformed element defined in the current element coordinate system may always be much smaller than unity. Thus the higher-order terms of deformation parameters in the element internal nodal forces may be neglected. However, in order to include the nonlinear coupling among the bending and stretching deformations, the terms up to the second order of deformation parameters and their spatial derivatives are retained in element deformation nodal forces by consistent second-order linearization ofδε11σ11 in (2.22). Here, only infinitesimal-free vibration is considered, thus only the terms up to the first order of time derivatives of deformation parameters and their spatial derivatives are retained in element inertia nodal forces by consistent first-order linearization of δrta in (2.22).From (2.6) and (2.11), the variation of ε11 in (2.13) may be expressed as(2.31)δε11=δεc-yδv,xx,δεc=1L(δuatGa+δubtGb),δv,xx=δubtNb′′.From (2.3), (2.6), and (2.12), δr the variation of r in (2.3) may be expressed as(2.32)δr={δr1,δr2,δr3}={δxp-yδv,x,δv-yv,xδv,x,0},δxp=δuatNa+xLδubtGb-∫0xv,xδv,xdx,δv,x=δubtNb′.Substituting (2.15)–(2.21) and (2.31)–(2.32) into (2.22), using ∫ydA=0, neglecting the higher order terms, we may obtain(2.33)faD=EAεcGa,(2.34)fbD=EI∫Nb′′v,xxdx+f12D∫Nb′v,xdx,(2.35)faI=ρA∫NaNatdxüa+Ω2ρAao1∫Nadx̲-Ω2ρA∫Na(Natua+x)dx̲-2ΩρA∫Nav̇dx,(2.36)fbI=ρA∫Nbv̈dx+ρI∫Nb′v̈,xdx+Ω2ρAao2∫Nbdx̲-Ω2ρA∫Nbvdx̲-Ω2ρI∫Nb′v′dx̲+2ΩρA∫NbNatdxu̇a.
where the range of integration for the integral ∫()dx in (2.34)–(2.36) is from 0 to L, A is the cross section area, I is moment of inertia of the cross section, aoi(i=1,2) are the xi components of ao in (2.17).The underlined terms in (2.35) and (2.36) are the inertia nodal force corresponding to the steady-state deformation induced by the constant rotation.


### 2.6. Element Matrices

The element matrices considered are element tangent stiffness matrix,mass matrix, centripetal stiffness matrix, and gyroscopic matrix. The element matrices may be obtained by differentiating the element nodal force vectors in (2.33)–(2.36) with respect to nodal parameters and time derivatives of nodal parameters.Using the direct stiffness method, the element tangent stiffness matrix may be assembled by the following submatrices:(2.37)kaa=∂faD∂ua=EALGaGat,kab=kbat=∂faD∂ub=0,kbb=∂fbD∂ub=EI∫Nb′′Nb′′tdx+f12D∫Nb′Nb′tdx.The element mass matrix may be assembled by the following submatrices:(2.38)maa=∂faI∂üa=ρA∫NaNatdx,mab=mbat=∂faI∂üb=0,mbb=∂fbI∂üb=ρA∫NbNbtdx+ρI(1-εc)2∫Nb′Nb′tdx.The element centripetal stiffness matrix may be assembled by the following submatrices:(2.39)kΩaa=∂faIΩ2∂ua=-ρA∫NaNatdx,kΩab=kΩabt=∂faIΩ2∂ub=0,kΩbb=∂fbIΩ2∂ub=-ρA∫NbNbtdx.The element gyroscopic matrix may be assembled by the following submatrices:(2.40)caa=∂faIΩ∂u̇a=0,cab=-cbat=∂faIΩ∂u̇b=-2ρA∫NaNbtdx,cbb=∂fbIΩ∂u̇b=0.


### 2.7. Equations of Motion

For convenience, the dimensionless variables defined in Table1 are used here.Table 1
Dimensionless variables.
VariablesDimensionless variablesCoordinatesx, Xo, Yox̅=xLT,X̅o=XoLT,Y̅o=Yo/LT,Timetτ=tLTEρLength of beam elementLL̅=L/LTArea moment of inertiaII̅=IALT2Radius of hubRR̅=R/LTDisplacementsu,vu̅=u/LT,v̅=v/LTspatial derivatives of displacementu′,u′′,v′,v′′u¯′=∂u¯∂x¯=u′,u¯′′=∂2u¯∂x¯2=LTu′′,v¯′=∂v¯∂x¯=v′,v¯′′=∂2v¯∂x¯2=LTv′′Time derivatives of displacementu̇,ü,v̇,v̈u¯̇=∂u¯∂τ=u̇ρE,u¯̈=∂2u¯∂τ2=LTρEü,v¯̇=∂v¯∂τ=v̇ρE,v¯̈=∂2v¯∂τ2=LTρEv̈Force and momentfij,mjf̅ij=fijEA,(i=1,2;j=1,2)Angular velocityΩk=ΩLTρ/ENatural frequencyωK=ωLTρ/EThe dimensionless nonlinear equations of motion for a rotating beam with constant angular velocity may be expressed by(2.41)φ=FD(Q⌢)+FI(k2,Q⌢,Q⌢̇,Q⌢̈)=0,(2.42)Q⌢=Qs+Q(τ),
where k and τ are dimensionless time and dimensionless angular speed of rotating beam, respectively, defined in Table 1. φ, FD, and FI are the dimensionless unbalanced force vector, the dimensionless deformation nodal force vector, and the dimensionless inertia nodal force vector of the structural system, respectively. FI and FD are assembled from the dimensionless element nodal force vectors, which are calculated using (2.33)–(2.36) and the dimensionless variables defined in Table 1 first in the current element coordinates and then transformed from element coordinate system to global coordinate system before assemblage using standard procedure. Q⌢ is the dimensionless nodal displacement vector of the rotating beam, Q⌢̇=∂Q⌢/∂τ and Q⌢̈=∂2Q⌢/∂τ2 are the dimensionless nodal velocity vector and the dimensionless nodal acceleration vector of the rotating beam, respectively, Qs is the dimensionless steady-state nodal displacement vector induced by constant dimensionless rotation speed k, and Q(τ) is the time-dependent dimensionless nodal displacements vector caused by the free vibration of the rotating beam. Here only infinitesimal vibration is considered.


### 2.8. Governing Equations for Steady-State Deformation

For the steady-state deformations,Q(τ)=0. Thus (2.41) can be reduced to nonlinear dimensionless steady-state equilibrium equations and expressed by(2.43)φ=FsD(Qs)+k2FsI(Qs)=0,
where FsD(Qs) and k2FsI(Qs) are the dimensionless deformation nodal force vector and the dimensionless inertia nodal force (the centrifugal force) vector of the structural system corresponding to the dimensionless steady-state nodal displacement vector Qs, respectively. k2FsI(Qs) is corresponding to the underlined terms of (2.35) and (2.36). Note that k2FsI(Qs) is deformation dependent.  Thus k2FsI(Qs) should be updated at each new configuration.Here, an incremental-iterative method based on the Newton-Raphson method is employed for the solution of nonlinear dimensionless steady-state equilibrium equations at different dimensionless rotation speedk. In this paper, a weighted Euclidean norm of the unbalanced force is employed for the equilibrium iterations and is given by(2.44)‖φ‖k2N‖FsI‖≤etol,
where N is number of the equations of the system and etol is a prescribed value of error tolerance. Unless otherwise stated, the error tolerance etol is set to 10-5 in this study.


### 2.9. Governing Equations for Free Vibration Measured from the Position of Steady-state Deformation

Substituting (2.42) into (2.41) and setting the first-order Taylor series expansion of the unbalanced force vector φ around Qs to zero, one may obtain the dimensionless governing equations for linear free vibration of the rotating beam measured from the position of the steady-state deformation as follows.(2.45)MQ̈+CQ̇+(K+k2KΩ)Q=0,
where M, C,K, and KΩ are dimensionless mass matrix, gyroscopic matrix, tangent stiffness matrix, and centripetal stiffness matrix of the rotating beam, respectively. M, C, K, and KΩ are assembled from the dimensionless element mass matrix, gyroscopic matrix, tangent stiffness matrix, and centripetal stiffness matrix, which are calculated using (2.37)–(2.40) and the dimensionless variables defined in Table 1 first in the current element coordinates and then transformed from element coordinate system to global coordinate system before assemblage using standard procedure.We will seek a solution of (2.45) in the form(2.46)Q=(QR+iQI)eiKτ,
wherei=-1, K and τ are dimensionless natural frequency of rotating beam and dimensionless time defined in Table 1, and QR and QI are real part and imaginary part of the vibration mode.Substituting (2.46) into (2.45), one may obtain a set of homogeneous equations expressed by(2.47)HZ=0,(2.48)H=H(K,k)=[K+k2KΩ-K2MkKCtkKCK+k2KΩ-K2M],(2.49)Z={QR,QI},
where H(K,k) denotes H being a function of Kand k. Note that H is a symmetric matrix.Equation (2.47) is a quadratic eigenvalue problem. For a nontrivial Z, the determinant of matrix H in (2.47) must be equal to zero. The values of K which make the determinant vanishes are called eigenvalues of matrix H. The bisection method is used here to find the eigenvalues. Note that when k=0, (2.47) will degenerate to a generalized eigenvalue problem.


## 2.1. Description of Problem

Consider an inclined uniform Euler beam of lengthLT rigidly mounted with an inclination angle α on the periphery of rigid hub with radius R rotating about its axis fixed in space at a constant angular speed Ω as shown in Figure 1. The axis of the rotating hub is perpendicular to one of the principal directions of the cross section of the beam. The deformation displacements of the beam are defined in an inertia rectangular Cartesian coordinate system which is coincident with a rotating rectangular Cartesian coordinate system rigidly tied to the hub.Here only axial and lagwise bending vibrations are considered. It is well known that the beam sustains a steady-state deformations (time-independent deformation displacements) induced by constant rotation [26]. In this study, large displacement and rotation with small strain are considered in the steady-state deformation. The vibration (time-dependent deformation displacements) of the beam is measured from the position of the steady-state deformation, and only infinitesimal-free vibration is considered. Note that the axial and lagwise vibrations, which are coupled due to the Coriolis effects and the lateral steady-state deformation, cannot be analyzed independently. Here the engineering strain and stress are used for the measure of the strain and stress.


## 2.2. Basic Assumptions

The following assumptions are made in derivation of the beam element behavior.(1)
The beam is prismatic and slender, and the Euler-Bernoulli hypothesis is valid.(2)
The unit extension of the centroid axis of the beam element is uniform.(3)
The deformation displacements and rotations of the beam element are small.(4)
The strains of the beam element are small.In conjunction with the corotational formulation and rotating frame method, the third assumption can always be satisfied if the element size is properly chosen. Thus, only the terms up to the second order of deformation parameters and their spatial derivatives are retained in element position vector, strain, and deformation nodal forces by consistent second-order linearization in this study.


## 2.3. Coordinate Systems

In order to describe the system, we define three sets of right-handed rectangular Cartesian coordinate systems.(1)
A rotating global set of coordinates,Xi(i=1,2,3) (see Figures 1 and 2); the coordinates rotate about the hub axis at a constant angular speed Ω as shown in Figure 1. The origin of this coordinate system is chosen to be the intersection of the centroid axes of the hub and the undeformed beam. The X1 axis is chosen to coincide with the centroid axis of the undeformed beam, and the X2 and X3 axes are chosen to be the principal directions of the cross section of the beam at the undeformed state. The direction of the axis of the rotating hub is parallel to the X3 axis. The nodal coordinates, nodal deformation displacements, absolute nodal velocity, absolute nodal acceleration, and equations of motion of the system are defined in terms of an inertia global coordinate system which is coincident with the rotating global coordinate system.(2)
Element coordinates;xi(i=1,2,3) (see Figure 2), a set of element coordinates is associated with each element, which is constructed at the current configuration of the beam element. The coordinates rotate about the hub axis at a constant angular speed Ω. The origin of this coordinate system is located at the element node 1, the centroid of the end section. The  x1 axis is chosen to pass through two end nodes of the element; the directions of the x2 and x3 axes are chosen to coincide with the principal direction of the cross section in the undeformed state. Because only the displacements in X1X2 plane are considered, the directions of x3 axis and X3 axis are coincident. The position vector, deformations, absolute velocity, absolute acceleration, internal nodal forces, stiffness matrices, and inertia matrices of the elements are defined in terms of an inertia element coordinate system which is coincident with the rotating element coordinate system.Figure 2
Coordinate systems.In this study, the direction of the axis of the rotating hub is parallel to theX3 axis and only the displacements in X1X2 plane are considered. Thus, the angular velocity of the hub referred to the global coordinates may be given by(2.1)ΩG={00Ω},
where the symbol {} denotes a column matrix, which is used through the paper.


## 2.4. Kinematics of Beam Element

LetQ (Figure 3) be an arbitrary point in the beam element and P the point corresponding to Q on the centroid axis. The position vector of point Q in the undeformed configurations referred to the current element coordinate system may be expressed as(2.2)r0={x,y,z}.Figure 3
Kinematics of Euler beam.Using the approximationcosθ≈1-(1/2)θ2, sinθ≈θ, and (1+εc)≈1, retaining all terms up to the second order, the position vector of point Q in the deformed configurations referred to the current element coordinate system may be expressed as(2.3)r={r1,r2,r3}={xp-yθ,y(1-12θ2)+v,z},(2.4)θ≈sin⁡θ=∂v(x,t)∂s=∂v(x,t)∂x∂x∂s=v′1+εc≈v′,(2.5)εc=∂s∂x-1,
where xp(x,t) and v(x,t) are the x1 and x2 coordinates of point P, respectively, in the deformed configuration, t is time, θ=θ(x,t) is the angle counterclockwise measured from x1 axis to the tangent of the centroid axis of the deformed beam, εc is the unit extension of the centroid axis, and s is the arc length of the deformed centroid axis measured from node 1 to point P. In this paper, ()′ denotes (),x=∂()/∂x.Here, the lateral deflection of the centroid axis,v(x,t), is assumed to be the Hermitian polynomials of x and may be expressed by(2.6)v(x,t)={N1,N2,N3,N4}t{v1,v1′,v2,v2′}=Nbtub,
where vj=vj(t) and vj′=vj′(t)(j=1,2) are nodal values of v and v,x, respectively, at nodes j. Note that, due to the definition of the element coordinates, the values of vj(j=1,2) are zero. However, their variations and time derivatives are not zero. Ni(i=1–4) are shape functions and are given by(2.7)N1=14(1-ξ)2(2+ξ),N2=L8(1-ξ2)(1-ξ),N3=14(1+ξ)2(2-ξ),N4=L8(-1+ξ2)(1+ξ),(2.8)ξ=-1+2xL,
where L is the length of the undeformed beam element.Making use of assumptionsv,x≪1 and εc≪1, the relationship between xp(x,t), v(x,t), and x in (2.3) may be approximated by(2.9)xp(x,t)=u1+∫0x(1+εc-12v,x2)dx,
where u1 is the displacement of node 1 in the x1 direction. Note that due to the definition of the element coordinate system, the value of u1 is equal to zero. However, the variation and time derivatives of u1 are not zero.The axial displacements of the centroid axis may be determined from the lateral deflections and the unit extension of the centroid axis using (2.9).From (2.9), one may obtain
(2.10)l=L+u2-u1=xc(L,t)-xc(0,t)=∫0L(1+εc-12v,x2)dx
in which ℓ is the current chord length of the centroid axis of the beam element and u2 is the displacement of node 2 in the x1 direction. Using the assumption of uniform extension of the centroid axis and (2.10), εc in (2.10) maybe expressed by(2.11)εc=1L(Gatua+12Gbtub),Ga={-1,1},ua={u1,u2},Gb=∫0LNb′v,xdx.Substituting (2.11) into (2.9), one may obtain(2.12)xp(x,t)=Natua+x+x2LGbtub-12∫0xv,x2dx,Na={1-ξ2,1+ξ2}.From (2.3) and the definition of engineering strain [32, 33], making use of the assumption of small strain, and retaining the terms up to the second order of deformation parameters, the engineering strain in the Euler beam may be approximated by(2.13)ε11=εc-yv′′.The absolute velocity and acceleration vectors of pointQ in the beam element may be expressed as(2.14)v={v1,v2,v3}=vo+Ω×r+ṙ,(2.15)a={a1,a2,a3}=ao+Ω̇×r+Ω×(Ω×r)+2Ω×ṙ+r̈,(2.16)vo=Ω×rAo,(2.17)ao={ao1,ao2,ao3}=Ω×(Ω×rAo),(2.18)Ω=AGEtΩG,(2.19)rAo=AGEtrAoG,(2.20)rAoG=rAO+rOoG={Rcos⁡α+Xo,-Rsinα+Yo,0},
where r is the position vector of point Q given in (2.3) referred to the current moving element coordinate system, the symbol ()̇ denotes time derivative, Ω is the vector of angular velocity referred to the current inertia element coordinates, ΩG is the angular velocity of the hub referred to the global coordinates given in (2.1), AGE is the transformation matrix between the current global coordinates and the current element coordinates, vo and ao are the absolute velocity and absolute acceleration of point o, the origin of the current element coordinates, Xo and Yo are coordinates of point o referred to the current global coordinates, R is the radius of the hub, and α is inclination angle of the rotating beam.  Ω×(Ω×r)and 2Ω×ṙ are centripetal acceleration and Coriolis acceleration, respectively. ṙ and r̈ are the velocity and acceleration of point Q relative to the current moving element coordinates. From (2.3), (2.11) and (2.12), ṙ and r̈ may be expressed as
(2.21)ṙ={ṙ1,ṙ2,ṙ3}={ẋp-yv̇,x,v̇-yv̇,xv,x,0},r̈={r̈1,r̈2,r̈3}={ẍp-yv̈,x,v̈-yv̇,x2-yv̈,xv,x,0},ẋp=Natu̇a+xLGbtu̇b-∫0xv,xv̇,xdx,ε̇c=1L(Gatu̇a+Gbtu̇b),ẍp=Natüa+xL(Gbtüb+Ġbtu̇b)-∫0x(v,xv̈,x+v̇,xv̇,x)dx,ε̈c=1L(Gatüa+Ġbtu̇b+Gbtüb).Note that the current element coordinates constructed at the current configuration of the beam element rotate about the hub axis at the angular velocity of the hub. Thus, the centripetal acceleration and Coriolis acceleration corresponding to the inertia forces of the rotating beam are unique. For nonrotating beam,Ω=0 and ṙ and r̈ are the absolute velocity and acceleration referred to the current element coordinate.


## 2.5. Element Nodal Force Vector

Letuj,δvj, and δvj′(j=1,2) denote the virtual displacements in the x1 and x2 directions of the current inertia element coordinates, and virtual rotations applied at the element nodes j. The element nodal force corresponding to virtual nodal displacements δuj, δvj, and δvj′(j=1,2) are fij, the forces in the xi(i=1,2) directions, and mj moments about the x3 axis, at element local nodes j.The element nodal force vector is obtained from the d’Alembert principle and the virtual work principle in the current inertia element coordinates. The virtual work principle requires that(2.22)δuatfa+δubtfb=∫V(δε11σ11+ρδrta)dV,(2.23)δua={δu1,δu2},(2.24)δub={δv1,δv1′,δv2,δv2′},(2.25)fa=faD+faI={f11,f12},(2.26)fb=fbD+fbI={f21,m1,f22,m2},(2.27)faD={f11D,f12D},(2.28)fbD={f21D,m1D,f22D,m2D},(2.29)faI,={f11I,f22I},(2.30)fbI,={f21I,m1I,f22I,m2I},
where fi(i=a,b) are the generalized force vectors corresponding to δua and δub, respectively, fiD and fiI(i=a,b) are element deformation nodal force vector and inertia nodal force vector corresponding to fi, respectively, V is the volume of the undeformed beam element, and δε11 is the variation of ε11 in (2.13) corresponding to δua and δub. σ11 is the engineering stress. For linear elastic material, σ11=Eε11, where E is Young’s modulus. ρ is the density, δr is the variation of r in (2.3) (referred to the current inertia element coordinate system) corresponding to δua and δub, and a is the absolute acceleration in (2.15).If the element size is chosen to be sufficiently small, the values of the deformation parameters of the deformed element defined in the current element coordinate system may always be much smaller than unity. Thus the higher-order terms of deformation parameters in the element internal nodal forces may be neglected. However, in order to include the nonlinear coupling among the bending and stretching deformations, the terms up to the second order of deformation parameters and their spatial derivatives are retained in element deformation nodal forces by consistent second-order linearization ofδε11σ11 in (2.22). Here, only infinitesimal-free vibration is considered, thus only the terms up to the first order of time derivatives of deformation parameters and their spatial derivatives are retained in element inertia nodal forces by consistent first-order linearization of δrta in (2.22).From (2.6) and (2.11), the variation of ε11 in (2.13) may be expressed as(2.31)δε11=δεc-yδv,xx,δεc=1L(δuatGa+δubtGb),δv,xx=δubtNb′′.From (2.3), (2.6), and (2.12), δr the variation of r in (2.3) may be expressed as(2.32)δr={δr1,δr2,δr3}={δxp-yδv,x,δv-yv,xδv,x,0},δxp=δuatNa+xLδubtGb-∫0xv,xδv,xdx,δv,x=δubtNb′.Substituting (2.15)–(2.21) and (2.31)–(2.32) into (2.22), using ∫ydA=0, neglecting the higher order terms, we may obtain(2.33)faD=EAεcGa,(2.34)fbD=EI∫Nb′′v,xxdx+f12D∫Nb′v,xdx,(2.35)faI=ρA∫NaNatdxüa+Ω2ρAao1∫Nadx̲-Ω2ρA∫Na(Natua+x)dx̲-2ΩρA∫Nav̇dx,(2.36)fbI=ρA∫Nbv̈dx+ρI∫Nb′v̈,xdx+Ω2ρAao2∫Nbdx̲-Ω2ρA∫Nbvdx̲-Ω2ρI∫Nb′v′dx̲+2ΩρA∫NbNatdxu̇a.
where the range of integration for the integral ∫()dx in (2.34)–(2.36) is from 0 to L, A is the cross section area, I is moment of inertia of the cross section, aoi(i=1,2) are the xi components of ao in (2.17).The underlined terms in (2.35) and (2.36) are the inertia nodal force corresponding to the steady-state deformation induced by the constant rotation.


## 2.6. Element Matrices

The element matrices considered are element tangent stiffness matrix,mass matrix, centripetal stiffness matrix, and gyroscopic matrix. The element matrices may be obtained by differentiating the element nodal force vectors in (2.33)–(2.36) with respect to nodal parameters and time derivatives of nodal parameters.Using the direct stiffness method, the element tangent stiffness matrix may be assembled by the following submatrices:(2.37)kaa=∂faD∂ua=EALGaGat,kab=kbat=∂faD∂ub=0,kbb=∂fbD∂ub=EI∫Nb′′Nb′′tdx+f12D∫Nb′Nb′tdx.The element mass matrix may be assembled by the following submatrices:(2.38)maa=∂faI∂üa=ρA∫NaNatdx,mab=mbat=∂faI∂üb=0,mbb=∂fbI∂üb=ρA∫NbNbtdx+ρI(1-εc)2∫Nb′Nb′tdx.The element centripetal stiffness matrix may be assembled by the following submatrices:(2.39)kΩaa=∂faIΩ2∂ua=-ρA∫NaNatdx,kΩab=kΩabt=∂faIΩ2∂ub=0,kΩbb=∂fbIΩ2∂ub=-ρA∫NbNbtdx.The element gyroscopic matrix may be assembled by the following submatrices:(2.40)caa=∂faIΩ∂u̇a=0,cab=-cbat=∂faIΩ∂u̇b=-2ρA∫NaNbtdx,cbb=∂fbIΩ∂u̇b=0.


## 2.7. Equations of Motion

For convenience, the dimensionless variables defined in Table1 are used here.Table 1
Dimensionless variables.
VariablesDimensionless variablesCoordinatesx, Xo, Yox̅=xLT,X̅o=XoLT,Y̅o=Yo/LT,Timetτ=tLTEρLength of beam elementLL̅=L/LTArea moment of inertiaII̅=IALT2Radius of hubRR̅=R/LTDisplacementsu,vu̅=u/LT,v̅=v/LTspatial derivatives of displacementu′,u′′,v′,v′′u¯′=∂u¯∂x¯=u′,u¯′′=∂2u¯∂x¯2=LTu′′,v¯′=∂v¯∂x¯=v′,v¯′′=∂2v¯∂x¯2=LTv′′Time derivatives of displacementu̇,ü,v̇,v̈u¯̇=∂u¯∂τ=u̇ρE,u¯̈=∂2u¯∂τ2=LTρEü,v¯̇=∂v¯∂τ=v̇ρE,v¯̈=∂2v¯∂τ2=LTρEv̈Force and momentfij,mjf̅ij=fijEA,(i=1,2;j=1,2)Angular velocityΩk=ΩLTρ/ENatural frequencyωK=ωLTρ/EThe dimensionless nonlinear equations of motion for a rotating beam with constant angular velocity may be expressed by(2.41)φ=FD(Q⌢)+FI(k2,Q⌢,Q⌢̇,Q⌢̈)=0,(2.42)Q⌢=Qs+Q(τ),
where k and τ are dimensionless time and dimensionless angular speed of rotating beam, respectively, defined in Table 1. φ, FD, and FI are the dimensionless unbalanced force vector, the dimensionless deformation nodal force vector, and the dimensionless inertia nodal force vector of the structural system, respectively. FI and FD are assembled from the dimensionless element nodal force vectors, which are calculated using (2.33)–(2.36) and the dimensionless variables defined in Table 1 first in the current element coordinates and then transformed from element coordinate system to global coordinate system before assemblage using standard procedure. Q⌢ is the dimensionless nodal displacement vector of the rotating beam, Q⌢̇=∂Q⌢/∂τ and Q⌢̈=∂2Q⌢/∂τ2 are the dimensionless nodal velocity vector and the dimensionless nodal acceleration vector of the rotating beam, respectively, Qs is the dimensionless steady-state nodal displacement vector induced by constant dimensionless rotation speed k, and Q(τ) is the time-dependent dimensionless nodal displacements vector caused by the free vibration of the rotating beam. Here only infinitesimal vibration is considered.


## 2.8. Governing Equations for Steady-State Deformation

For the steady-state deformations,Q(τ)=0. Thus (2.41) can be reduced to nonlinear dimensionless steady-state equilibrium equations and expressed by(2.43)φ=FsD(Qs)+k2FsI(Qs)=0,
where FsD(Qs) and k2FsI(Qs) are the dimensionless deformation nodal force vector and the dimensionless inertia nodal force (the centrifugal force) vector of the structural system corresponding to the dimensionless steady-state nodal displacement vector Qs, respectively. k2FsI(Qs) is corresponding to the underlined terms of (2.35) and (2.36). Note that k2FsI(Qs) is deformation dependent.  Thus k2FsI(Qs) should be updated at each new configuration.Here, an incremental-iterative method based on the Newton-Raphson method is employed for the solution of nonlinear dimensionless steady-state equilibrium equations at different dimensionless rotation speedk. In this paper, a weighted Euclidean norm of the unbalanced force is employed for the equilibrium iterations and is given by(2.44)‖φ‖k2N‖FsI‖≤etol,
where N is number of the equations of the system and etol is a prescribed value of error tolerance. Unless otherwise stated, the error tolerance etol is set to 10-5 in this study.


## 2.9. Governing Equations for Free Vibration Measured from the Position of Steady-state Deformation

Substituting (2.42) into (2.41) and setting the first-order Taylor series expansion of the unbalanced force vector φ around Qs to zero, one may obtain the dimensionless governing equations for linear free vibration of the rotating beam measured from the position of the steady-state deformation as follows.(2.45)MQ̈+CQ̇+(K+k2KΩ)Q=0,
where M, C,K, and KΩ are dimensionless mass matrix, gyroscopic matrix, tangent stiffness matrix, and centripetal stiffness matrix of the rotating beam, respectively. M, C, K, and KΩ are assembled from the dimensionless element mass matrix, gyroscopic matrix, tangent stiffness matrix, and centripetal stiffness matrix, which are calculated using (2.37)–(2.40) and the dimensionless variables defined in Table 1 first in the current element coordinates and then transformed from element coordinate system to global coordinate system before assemblage using standard procedure.We will seek a solution of (2.45) in the form(2.46)Q=(QR+iQI)eiKτ,
wherei=-1, K and τ are dimensionless natural frequency of rotating beam and dimensionless time defined in Table 1, and QR and QI are real part and imaginary part of the vibration mode.Substituting (2.46) into (2.45), one may obtain a set of homogeneous equations expressed by(2.47)HZ=0,(2.48)H=H(K,k)=[K+k2KΩ-K2MkKCtkKCK+k2KΩ-K2M],(2.49)Z={QR,QI},
where H(K,k) denotes H being a function of Kand k. Note that H is a symmetric matrix.Equation (2.47) is a quadratic eigenvalue problem. For a nontrivial Z, the determinant of matrix H in (2.47) must be equal to zero. The values of K which make the determinant vanishes are called eigenvalues of matrix H. The bisection method is used here to find the eigenvalues. Note that when k=0, (2.47) will degenerate to a generalized eigenvalue problem.


## 3. Numerical Examples

To verify the accuracy of the present method and to investigate the steady deformation and the natural frequencies of rotating-inclined beams with different inclination angleα, dimensionless radius of the hub R̅, and slenderness ratios η=LTA/I at different dimensionless angular velocities k, several dimensionless numerical examples are studied here.For simplicity, only the uniform beam with rectangular cross section is considered here. The maximum steady-state axial strainεmax of rotating beam is the sum of the maximum steady-state membrane strain εcmax and bending strain εbmax, which occur at the root of the rotating beam. In practice, rotating structures are designed to operate in the elastic range of the materials. Thus, it is considered that εmax≤εy (say 0.01) in this study. At the same dimensionless angular speed k, εmax are different for rotating beams with different η, α, and R̅.  Thus, the allowable k are different for rotating beams with different η, α, and R̅ in this study.To investigate the effect of the lateral deflection on the steady-state deformation and the natural frequency of rotating Euler beams, here cases with and without considering the lateral deflection are considered. The corresponding elements are referred to as EA element and EB element, respectively. For EA element, all terms in (2.33)–(2.40) are considered; for EB element, all terms in (2.33)–(2.40) are considered except the underlined terms in (2.36), which are the lateral inertia nodal force corresponding to the steady-state deformation induced by the constant rotation. In this section, vtip/LT denotes the dimensionless lateral tip deflection of the steady-state deformation; Ki denotes the ith dimensionless natural frequency of the rotating beam and denote that the corresponding vibration mode is lateral vibration at k=0; in all tables, the entries with “(a)” denotes that the corresponding vibration mode is axial vibration at k=0.The example first considered is the rotating-inclined beams with dimensionless radius of the hubR̅=1.5, inclination angle α=0°,5°,30°,90°, and slenderness ratios η=20,1000. The present results are shown in Tables 2 and 3 together with some results available in the literature. In Tables 2 and 3,  EAn, n=10,50,100, denote that n equal EA elements are used for discretization, and LAS denotes the linear analytical solution of the steady-state deformation. It can be seen that for higher natural frequencies of lateral vibration, the discrepancy between the present results and the analytical solutions given in [34], in which the rotary inertia is not considered, increases with decrease of the slenderness ratio. It seems that the effect of the rotary inertia on the higher natural frequencies of the Euler beam is not negligible when the slenderness ratio is small. It can be seen from Tables 2 and 3 that the differences between the results of EA50 and EA 100 are negligible for all cases studied. Thus, in the rest of the section, all numerical results are obtained using 50 equal elements. For α=0, and k≠0, the steady-state deformation is axial deformation only as expected. The analytical solution of the maximum steady-state membrane strain εcmax=k2(R̅cosα+1/2) given in [15] and the linear solution are identical. It can be seen that at the same dimensionless angular speed k, εcmax is independent of the slenderness ratio η.  Thus, for α=0, the allowable k is limited by εcmax and is the same for the rotating beam with different slenderness ratio η. Very good agreement is observed between the natural frequencies obtained by the present study and those given in [24], which are obtained using the power series method. It can be seen from Table 3 that for slenderness ratio η = 1000, with increase of the inclination angle α, the values of εbmax and vtip/LT increase significantly and the value of the allowable dimensionless angular speed k decreases significantly. Comparing εbmax and vtip/LT of EA with the results of linear analytical solution, respectively, it is found that the difference between the results of EA and LAS is insignificant for η=20 but is remarked for η=1000.  These may be explained as follows. The centrifugal stiffening effect is significant for slender beam, and the lateral component of the centrifugal force in the rotating inclined beam decreases with the increase of the steady-state lateral deflection.Table 2
Comparison of results for different cases (η=20, R̅=1.5).
αkεcmax(10-3)εbmax(10-3)vtip/LT(10-3)K1K2K3(a)K4K5(a)K6K7(a)EA10000.1747881.059571.572412.824954.756105.195468.00214EA50000.1747871.059531.570862.824314.714135.191207.862060EA100000.1747871.059531.570812.824314.712835.191197.85600[24]000.174791.059531.570802.824314.712395.19119—0°[34]000.175801.101721.570803.084864.712396.04510—EA106.9330900.1986161.087561.576152.853334.757295.223848.02928EA507.1549200.1985141.087261.576162.852434.715345.219317.862740.06EA1007.1821000.1985111.087261.574572.852424.714035.219307.85669[24]7.2000000.198621.087601.574552.852764.713605.21962—LAS7.2000000———————EA101.726801.930985.47630.1810491.066611.573352.832064.756395.202568.008895°0.03EA501.781951.935465.47699.1810211.066511.571802.831364.714435.198237.86221EA1001.788701.935605.47701.1810201.066511.571752.831364.713125.198227.85616LAS1.794862.037945.88301———————EA10.1732981.290083.72294.1754101.060281.572522.825674.756135.196198.0028130°0.01EA50.1786151.292243.72299.1754071.060241.570972.825034.714165.191917.86207EA100.1792641.292313.72300.1754071.060241.570922.825034.712855.191907.85601LAS.1799041.299043.75000———————EA10.05003452.593647.49504.1748361.059781.572532.825204.756125.195738.0022990°0.01EA50.05003842.597847.49506.1748351.059741.570982.824564.714155.191457.86205EA100.05002162.597977.49507.1748351.059741.570932.824564.712845.191447.85599LAS.05000002.598077.50000———————Table 3
Comparison of results for different cases (η=1000, R̅=1.5).
αkεcmax(10-3)εbmax(10-3)vtip/LTK1(10-2)K2(10-1)K3(10-1)K4K5K6K7EA10000.351601.220349.617105.121008.200340.300117.421052EA50000.351601.220341.616949.120893.199838.298509.4169030EA100000.351601.220341.616948.120893.199837.298506.416896[24]000.352.2203.6169.12089.19984.29851—0°[34]000.3516.22034.616972.120902———EA106.93309009.004572.501864.13423.591446.784725.9929271.21760EA507.15492008.962392.474244.06068.580524.771309.9761201.193650.06EA1007.18210008.961522.473124.05756.580088.770833.9756341.19316[24]7.20000008.9522.47084.0536.57955.77017.97486—LAS7.2000000———————EA101.731133.88303.08351714.547141.274482.17658.323098.443024.577262.7266985°0.03EA501.783966.00526.08381944.533481.262202.15028.319777.439167.572464.719965EA1001.789366.20203.08382184.533201.261792.14942.319678.439068.572368.719873LAS1.79486101.89714.70753———————EA10.1171748.73588.4296881.29068.405580.836390.143462.221484.319586.43911030°0.008EA50.1143419.36150.4299791.28848.404156.836065.143631.221433.318254.434653EA100.1134109.38784.4299861.28840.404108.836056.143643.221458.318289.434691LAS.11513841.56926.00000———————EA10006325878.11012.747138.561367.232168.566051.113317.190889.289637.40972290°0.003EA50003882248.15298.747250.560585.232182.566299.113204.190326.287895.405359EA100.003517408.15396.747254.560558.232181.566306.113202.190322.287886.405342LAS.0045000011.69131.68750———————To investigate the effect of the lateral deflection on the steady-state deformation and the natural frequency of rotating-inclined beams, the cases with and without considering the lateral deflection are studied forη=70, R̅=1, and k=5/70. The present results are shown in Table 4. The results transcribed from the figure given in [21], in which the steady-state lateral deflection and the rotary inertia are not considered, are also shown in Table 4 for comparison. It can be seen from Table 4 that except α=0, the values of εbmax are much larger than the yield strain for most engineering materials at k=5/70. Thus the results in Table 4 are only displayed for the purpose of comparisons between the results of EB and those given in [21]. There is a very good agreement between the natural frequencies obtained using the EB element and those given in [21]. Although the comparisons are beyond the yield point of most engineering materials, results of EA and EB show that the differences between the cases with and without considering the lateral deflection become apparent for the rotating-inclined beam with large inclination angle α at high-dimensionless angular speed. It can be seen from Table 4 that the difference between the natural frequencies of EA and EB is not significant for small α, but the first natural frequency of EB is much smaller than that of EA for large α.  The natural frequencies of EA slightly decrease with increase of α, but those of EB significantly decrease with increase of α for α≥50°. These may be partially attributed to the factthat the decrease of the centrifugal stiffening effect of the rotating-inclined beam caused by the increase of the inclination angle is alleviated by the increase of lateral deflection induced by the lateral centrifugal force.Table 4
Dimensionless frequencies for rotating beam with different inclination angle  (η=70, R̅=1,k=5/70).
αεcmax(10-3)εbmaxvtip/LTK1K2EAEBEAEAEAEB[21]EAEB[21]0°7.615827.6157900.105565.105427.105.411754.410792.41810°7.532757.53893.021843.119537.105513.104869.105.411356.410001.41720°7.285567.31066.043509.236923.105359.103195.103.410160.407642.41430°6.880326.93792.064820.350059.105101.100399.100.408162.403758.41040°6.327006.43205.085602.456935.104737.0964721.096.405357.398421.40550°5.639275.80840.105681.555686.104264.0913941.091.401742.391733.39860°4.834255.08596.124890.644628.103679.0851262.085.397314.383830.39070°3.932144.28663.143064.722301.102976.0775919.077.392079.374876.38180°2.955753.43472.160043.787503.102151.0686418.068.386048.365073.37190°1.930102.55611.175673.839324.101193.0579597.057.379245.354659.361To investigate the effect of angular speed on the steady-state deformation and natural frequency of rotating beams with different slenderness ratios and inclination angles, the following cases are considered: slenderness ratioη=39,50,100,1000, inclination angle α=0°,5°,30°,90°, and dimensionless radius of the rotating hub R̅=1. Tables 5, 6, 7, and 8 tabulate the maximum steady-state membrane strain and bending strain, dimensionless lateral tip deflection, and first six dimensionless natural frequencies for different η. It can be seen from Tables 5–8 that the values of vtip/LT increase significantly with the increase of the dimensionless angular velocities k and slenderness ratio η.  However, the values of vtip/LT are very small for η=39 and 50. Because the stiffening effect of the centrifugal force is significant for slender beam, as expected, it can be seen from Table 8 that the lower natural frequencies of lateral vibration increase remarked with increase of the dimensionless angular speed for η=1000.Table 5
Dimensionless frequencies for rotating beam with different inclination angle (η=39, R̅=1).
αkεcmax(10-4)εbmax(10-3)vtip/LT(10-3)K1(10-1)K2K3K4K5K6(a)0000.900168.5590571.543251.57086 (a)2.963964.71413.0101.4899900.909817.5602831.544521.57097 (a)2.965284.71415.0205.9606000.938126.5639451.548301.57130 (a)2.969264.714210°.03013.413800.983359.5699961.554551.57190 (a)2.975864.71432.04023.8528001.04313.5783591.563021.57296 (a)2.985094.71449.05037.2822001.11488.5889351.571131.57704 (a)2.996914.71472.06053.7078001.19620.6016041.57358 (a)1.589363.011294.71502.005.371545.073310.412011.902582.5593631.543571.57089 (a)2.964294.71414.0101.48623.2899771.62160.909786.5602801.544491.57100 (a)2.965284.714155°.0153.34417.6407073.55360.921661.5618051.545951.57126 (a)2.966934.71416.0205.945591.111506.09535.938019.5639321.547841.57177 (a)2.969254.71416.0259.290751.685479.11211.958617.5666541.550031.57268 (a)2.972224.71415.03013.38002.3447312.4630.983170.5699621.552331.57416 (a)2.975854.71413.002.054292.067510.379956.900509.5591021.543301.57087 (a)2.964014.71413.004.217174.2696141.51631.901534.5592371.543411.57090 (a)2.964154.7141330°.006.488657.6050413.39858.903243.5594611.543551.57103 (a)2.964404.71412.008.8687631.071706.00953.905635.5597751.543621.57132 (a)2.964744.71409.0101.357521.666739.32547.908710.5601771.543511.57191 (a)2.965184.71403.002.019999.135091.760426.900208.5590761.543271.57087 (a)2.963984.71413.004.080009.5403623.04080.900333.5591321.543241.57098 (a)2.964044.7141290°.006.1800721.215796.83830.900554.5592241.542951.57140 (a)2.964154.71407.008.3202532.1613212.1478.900889.5593491.542111.57246 (a)2.964314.71394.010.5006313.3767818.9613.901363.5595061.540411.57447 (a)2.964514.71365Table 6
Dimensionless frequencies for rotating beam with different inclination angle (η=50,R̅=1).
αkεcmax(10-4)εbmax(10-3)vtip/LT(10-3)K1(10-1)K2K3K4 (a)K5K60000.702550.4378591.215301.570862.351763.82646.0101.4899800.714917.4394411.216961.570962.353493.82822.0205.9605800.750712.4441541.221921.571252.358683.833490°.03013.413700.806600.4518981.230141.571732.367313.84227.04023.852700.878442.4625161.241551.572412.379323.85452.05037.282000.962316.4758131.256051.573282.394673.87021.06053.7076001.05500.4915631.273521.574352.413293.88930.005.371544.093761.674835.705653.4382541.215711.570892.352193.82689.0101.48623.3682792.62905.714878.4394371.216951.570962.353493.828215°.0153.34418.8048925.66932.729982.4414021.219001.571112.355653.83040.0205.945601.377179.52433.750592.4441361.221841.571342.358693.83347.0259.290742.0559413.8991.776240.4476261.225471.571662.362583.83740.03013.37992.8135718.5209.806412.4518541.229871.572082.367343.84221.002.054292.086522.624206.702990.4379171.215361.570872.351823.82652.004.217177.3451912.48734.704308.4380911.215541.570882.352023.8267130°.006.488671.7733495.56121.706507.4383811.215831.570932.352343.82704.008.8688031.366659.79986.709585.4387851.216211.571032.352803.82749.0101.357602.1192115.1409.713543.4393011.216661.571202.353403.82806.002.019999.1731931.24980.702604.4378831.215321.570872.351793.82648.004.080020.6927674.99671.702772.4379561.215391.570902.351883.8265790°.006.1801291.5586611.2327.703078.4380731.215451.571032.352043.82670.008.3204192.7706519.9428.703561.4382321.215441.571352.352283.82688.010.5009974.3281431.1020.704273.4384261.215271.572002.352653.82709Table 7
Dimensionless frequencies for rotating beam with different inclination angle (η=100, R̅=1).
αkεcmax(10-4)εbmax(10-3)vtip/LT(10-3)K1(10-1)K2K3K4 (a)K5K60.00000000.351520.219989.6146021.200471.570861.97619.0101.4899800.375696.223163.6179631.204021.570961.97983.0205.9605800.439855.232418.6279261.214581.571241.990720°.03013.413700.528670.247050.6441401.231941.571722.00872.04023.852700.630838.266138.6660821.255781.572402.03361.05037.282000.740095.288755.6931271.285681.573272.06512.06053.707600.853195.314093.7246141.321161.574332.10288.005.371547.1841392.62891.357708.220785.6154411.201361.570891.97710.0101.48624.6885859.52242.375636.223154.6179401.203971.571051.979855°.0153.344161.4067918.5127.403666.227045.6220621.208261.571401.98446.0205.945402.2383827.6550.439745.232381.6277731.214191.571951.99093.0259.290103.1197035.8627.481914.239066.6350401.221761.572611.99922.03013.37864.0196042.7784.528564.246988.6438161.230961.573292.00929.002.054294.1725962.48732.352403.220106.6147251.200601.570871.97632.004.217197.6833269.79956.355052.220456.6150811.200961.570971.9767530°.006.4887561.5118821.5060.359467.221033.6156361.201451.571361.97753.008.8690032.6271336.9471.365637.221832.6163441.201971.572271.97873.0101.357903.9906555.3127.373528.222846.6171631.202411.573951.98043.002.020005.3463834.99669.351634.220038.6146511.200521.570891.97625.004.0801031.3853219.9425.352038.220178.6147481.200531.571271.9765490°.006.1804783.1149744.6708.352915.220392.6147361.200141.572881.97728.008.3212265.5283178.7975.354553.220651.6143731.198781.577041.97891.010.5020648.60794121.586.357326.220924.6133681.195811.585231.98201Table 8
Dimensionless frequencies for rotating beam with different inclination angle (η=1000, R̅=1).
αkεcmax(10-5)εbmax(10-3)vtip/LT(10-2)K1(10-2)K2(10-1)K3(10-1)K4K5K60.00000000.351601.220341.616949.120893.199838.298509.01014.8998001.32581.432338.886010.151855.233376.333736.02059.6058002.52842.7661721.38410.216189.309541.4190560°.030134.137003.736481.114021.92511.289167.400408.525942.040238.527004.944951.465972.47903.364653.496007.640880.050372.820006.152671.819833.03912.441165.593428.759068.060537.076007.359072.174823.60301.518237.691707.878747.0053.71469.7651475.97092.741097.289190.695401.129392.208701.307598.01014.86001.675687.339381.32573.432307.885745.151766.233190.3334325°.01533.44312.579667.791731.92549.5959111.12500.182037.268266.371976.02059.47133.469778.020502.52834.7661601.38402.216153.309441.418839.02592.95134.341898.159573.13220.9393071.65211.252155.353965.470793.030133.8895.193168.253603.736411.114011.92507.289150.400356.525817.002.5424101.3144916.2578.436039.231506.625674.121676.200527.29913830°.0042.157093.4549931.0603.628044.264068.655502.124160.202509.300711.0064.842035.6188137.4726.851755.312455.712896.129772.207598.305211.0088.612687.7850240.65361.08530.369468.788930.138176.215896.313013.0005.012524.2164073.11010.352489.220632.617142.120911.199852.298522.0010.050204.86079412.1582.357434.221288.616084.120749.199640.298279.0015.1120681.8952825.7173.371925.221920.610342.119952.198654.29718190°.0020.1941233.2165340.7205.400095.222722.599365.118425.196788.295121.0025.2907754.6887954.0272.440751.224723.587411.116669.194633.292734.0030.3983556.2130164.2905.489485.228876.578807.115177.192740.290595.0035.5160567.7438571.7152.542488.235452.575106.114128.191302.288900Figures4–6 show the deformed configurations, axial displacements, and lateral displacements for the steady-state deformation of rotating beams with η=100, α=90°, and η=1000, α=5°,90° at different dimensionless angular speeds. In Figures 4–6, the X1 and X2 coordinates of the deformed configurations of rotating beam are present at the same scale, and X10 denotes the global Lagrangian coordinate of the beam axis. Very large displacement and rotation are observed in Figure 6.The steady-state deformation of rotating beam, (a) deformed configuration, (b) axial displacement, and (c) lateral displacement (η=100, R̅=1, α=90°).
(a)(b)(c)The steady-state deformation of rotating beam, (a) deformed configuration, (b) axial displacement, and (c) lateral displacement (η=1000, R̅=1, α=5°).
(a)(b)(c)The steady-state deformation of rotating beam, (a) deformed configuration, (b) axial displacement, (c) lateral displacement (η=1000, R̅=1, α=90°).
(a)(b)(c)Figures7–10 show the first six vibration modes for rotating beams with =39, α=0°,5°, and  η=1000, α=5°,90° at different dimensionless angular speeds. In Figures 7–10, U and V denote the X1 and X2 components of the vibration mode, respectively. The definitions of U and V are given by
(3.1)U=(UR2+UI2)1/2sign⁡(sinϕu),sinϕu=UI(UR2+UI2)1/2,-π≤ϕu≤π,V=(VR2+VI2)1/2sign⁡(sinϕv),sinϕv=VI(VR2+VI2)1/2,-π≤ϕv≤π,sign⁡(x)={1forx>0,-1forx<0,
where UR and VR, and UI and VI are the X1 and X2 components of QR and QI, real part and imaginary part of the vibration mode given in (2.47), respectively.  ϕu and ϕv are phase angles. Mode j(j=1–6) denotes the vibration mode dominated by the vibration mode corresponding to the jth natural frequency of the nonrotating beam. It can be seen from Figures 7–10, and Tables 5 and 8 that all vibration modes shown in Figures 7–10 are lateral vibration at k=0, except the fourth and the sixth vibration modes of η=39. Mode 3 is the third bending vibration mode, and mode 4 is the first axial vibration mode for η=39. Figure 11 shows the third and the fourth dimensionless natural frequencies for the rotating beam with η=39 at different dimensionless angular velocities. Because the third and the fourth natural frequencies are relatively close, frequency veering phenomenon [24] induced by the Coriolis force and the centrifugal force is observed in Figure 11. It can be seen from Figures 7 and 8 that the coupling of the axial and lateral vibration modes is very significant. Due to the effect of Coriolis force and the lateral steady-state deformation, the axial and lateral vibrations of rotating beam should be coupled.  However, from the numerical results of this study, it is found that the coupling is negligible for rotating beam with small slenderness ratio if the corresponding natural frequencies of the axial and lateral vibrations are not close. Due to the steady-state lateral deformation, it can be seen from Figures 9 and 10 that when k≠0, all vibration modes consist of the X1 and X2 components. The difference between the vibration modes of rotating beam at different k is very significant for η=1000.The first six vibration mode shapes of a rotating beam (η=39, R̅=1, α=0°).
(a)(b)(c)(d)(e)(f)The first six vibration mode shapes of a rotating beam (η=39,R̅=1,α=5°).
(a)(b)(c)(d)(e)(f)The first six vibration mode shapes of a rotating beam (η=1000, R̅=1, α=5°).
(a)(b)(c)(d)(e)(f)The first six vibration mode shapes of a rotating beam (η=1000, R̅=1, α=90°).
(a)(b)(c)(d)(e)(f)Figure 11
The third and fourth natural frequencies verse the dimensionless angular velocity (α=0°, r=1, η=39).


## 4. Conclusions

In this paper, a corotational finite element formulation combined with the rotating frame method and numerical procedure are proposed to derive the equations of motion for a rotating-inclined Euler beam at constant angular velocity. The element deformation and inertia nodal forces are systematically derived by the virtual work principle, the d’Alembert principle, and consistent second-order linearization of the fully geometrically nonlinear beam theory in the current element coordinates. The equations of motion of the system are defined in terms of an inertia global coordinate system which is coincident with a rotating global coordinate system rigidly tied to the rotating hub, while the total strains in the beam element are measured in an inertia element coordinate system which is coincident with a rotating element coordinate system constructed at the current configuration of the beam element. The rotating element coordinates rotate about the hub axis at the angular speed of the hub. The steady-state deformation and the natural frequency of infinitesimal-free vibration measured from the position of the corresponding steady-state deformation are investigated for rotating-inclined Euler beams with different inclination angles, slenderness ratios, and angular speeds of the hub.The results of dimensionless numerical examples demonstrate the accuracy and efficiency of the proposed method. The present results show that the geometrical nonlinearities that arise due to steady-state lateral and axial deformations should be considered for the natural frequencies of the inclined-rotating beams. Due to the effect of the centrifugal stiffening, the lower dimensionless natural frequencies of lateral vibration increase remarked with increase of the dimensionless angular speed for slender beam. The decrease of the centrifugal stiffening effect of the rotating inclined beam caused by the increase of the inclination angle is alleviated by the increase of lateral deflection induced by the lateral centrifugal force. Due to effect of the Coriolis force and centrifugal stiffening, frequency veering phenomenon is observed when inclination angleα=0°, and two natural frequencies corresponding to axial vibration and lateral vibration are close.Finally, it may be emphasized that, although the proposed method is only applied to the two dimensional rotating cantilever beams with inclination angle here, the present method can be easily extended to three dimensional rotating beams with precone and setting angle.



---

*Source: 146505-2011-08-04.xml*