# Free Vibration Analysis of a Rotating Composite Shaft Using thep-Version of the Finite Element Method

Authors: A. Boukhalfa; A. Hadjoui; S. M. Hamza Cherif
Journal: International Journal of Rotating Machinery
 (2008)

---

## Abstract

This paper is concerned with the dynamic behavior of the rotating composite shaft on rigid bearings. Ap-version, hierarchical finite element is employed to define the model. A theoretical study allows the establishment of the kinetic energy and the strain energy of the shaft, necessary to the result of the equations of motion. In this model the transverse shear deformation, rotary inertia and gyroscopic effects, as well as the coupling effect due to the lamination of composite layers have been incorporated. A hierarchical beam finite element with six degrees of freedom per node is developed and used to find the natural frequencies of a rotating composite shaft. A program is elaborate for the calculation of the eigenfrequencies and critical speeds of a rotating composite shaft. To verify the present model, the critical speeds of composite shaft systems are compared with those available in the literature. The efficiency and accuracy of the methods employed are discussed.


---

## Body

## 1. Introduction

The
application of composite shafts has come a long way from early low-speed
automotive driveshafts to helicopter tail rotors operating above the second
critical speed. With operation at supercritical speeds, a substantial amount of
payoffs and net system weight reductions are possible. At the same time, the
rotordynamic aspects assume more importance, and detailed analysis is required.
There are some technological problems associated with implementation, such as
joints with bearings, affixing of lumped masses, couplings, provision of
external damping, and so forth. The solutions proposed are just adequate, but
require substantial refinements, which might explain some of the differing
experiences of various authors.Zinberg
and Symonds [1] described a boron/epoxy composite tail rotor driveshaft for a
helicopter. The critical speeds were determined using equivalent modulus beam
theory (EMBT), assuming the shaft to be a thin-walled circular tube simply
supported at the ends. Shear deformation was not taken into account. The shaft
critical speed was determined by extrapolation of the unbalance response curve
which was obtained in the subcritical region.Dos
Reis et al. [2] published analytical investigations on thin-walled layered
composite cylindrical tubes. In part III of the series of publications, the
beam element was extended to formulate the problem of a rotor supported on
general eight coefficient bearings. Results were obtained for shaft
configuration of Zinberg and Symmonds. The authors have shown that
bending-stretching coupling and shear-normal coupling effects change with
stacking sequence, and alter the frequency values. Gupta and Singh [3] studied
the effect of shear-normal coupling on rotor natural frequencies and modal
damping. Kim and Bert [4] have formulated the problem of determination of
critical speeds of a composite shaft including the effects of bending-twisting
coupling. The shaft was modelled as a Bresse-Timoshenko beam. The shaft
gyroscopics have also been included. The results compare well with Zinberg's
rotor [1]. In another study, Bert and Kim [5] have analysed the dynamic
instability of a composite drive shaft subjected to fluctuating torque and/or
rotational speed by using various thin shell theories. The rotational effects
include centrifugal and Coriolis forces. Dynamic instability regions for a long
span simply supported shaft are presented.M.-Y. Chang et al. [6] published the vibration behaviours
of the rotating composite shafts. In  the model, the transverse shear
deformation, rotary inertia, and gyroscopic effects, as well as the coupling
effect due to the lamination of composite layers have been incorporated. The
model based on a first-order shear deformable beam theory (continuum-based
Timoshenko beam theory).C.-Y.
Chang et al. [7] published the vibration analysis of rotating composite shafts
containing randomly oriented reinforcements. The Mori-Tanaka mean-field theory
is adopted here to account for the interaction at the finite concentrations of
reinforcements in the composite material.Additional
recent work on composite shafts dealing with both the theoretical and
experimental aspects was reported by Singh [8], Gupta and Singh [3], and Singh
and Gupta [9, 10]. Rotordynamic formulation based on equivalent modulus beam theory
was developed for a composite rotor with a number of lumped masses and
supported on general eight coefficient bearings. A layerwise beam theory was derived
by Gupta and Singh  [3] from an available shell theory, with a layerwise
displacement field, and was then extended to solve a general composite
rotordynamic problem. The conventional rotor dynamic parameters as well as critical
speeds, natural frequencies, damping factors, unbalance response, and threshold
of stability were analysed in detail and results from the formulations based on
the two theories, namely, the equivalent modulus beam theory (EMBT) and
layerwise beam theory (LBT) were compared [9, 10]. The experimental rotordynamic
studies carried by Singh and Gupta [11, 12] were conducted on two filament wound
carbon/epoxy shafts with constant winding angles (±45° and ±60°). Progressive balancing had to be carried
out  to enable the shaft to traverse through the first critical speed. Inspite of the very different shaft configurations used, the authors' have
shown that bending-stretching coupling and shear-normal coupling effects change
with stacking sequence, and alter the frequency values.Some
practical aspects such as effect of shaft disc angular misalignment,
interaction between shaft bow, which is common in composite shafts and rotor
unbalance, and an unsuccessful operation of a composite rotor with an external
damper were discussed and reported by Singh and Gupta [11]. The Bode and
cascade plots were generated and orbital analysis at various operating speeds
was performed. The experimental critical speeds showed good correlation with
the theoretical prediction.This paper deals with thep-version,
hierarchical finite element method applied to free vibration analysis of rotating
composite shafts. The hierarchical concept for finite element
shape functions has been investigated during the past 25 years. Babuska et al. [13] established a theoretical basis for p-elements,
where the mesh keeps unchanged and the polynomial degree of the shape functions
is increased; however, in the standard h-version of the finite element
method the mesh is refined to achieve convergence and the polynomial degree of
the shape functions remains unchanged. Since then, standard forms of the
hierarchical shape functions have been represented in the literature elsewhere;
see for instance [14, 15].Meirovitch and Baruh [16] and Zhu [17] have shown that the hierarchical finite element
method yields a better accuracy than the h-version for eigenvalues
problems. The  hierarchical shape functions  used by Bardell [18] are based on integrated Legendre orthogonal polynomials; the symbolic
computing is used to calculate the mass and stiffness matrices of beams and
plates. Côté and Charron [19]
give the selection of p-version
shape functions for plate vibration analysis.In the presented composite
shaft model, the Timoshenko theory will be adopted. It is the purpose of the
present work to study dynamic characteristics such as natural frequencies,
whirling frequencies, and the critical speeds of the rotating composite shaft. In the model, the transverse shear
deformation, rotary inertia, and gyroscopic effects, as well as the coupling
effect due to the lamination of composite layers have been incorporated. To
determine the rotating shaft system's responses, the hierarchical finite
element method with trigonometric shape functions [20, 21] is used here to
approximate the governing equations by a system of ordinary differential equations.


## 2. Equations of Motion

## 2.1. Kinetic and Strain Energy Expressions

The
shaft is modelled as a Timoshenko beam, that is, first-order shear deformation
theory with rotary inertia and gyroscopic effect is used. The shaft rotates at
constant speed about its longitudinal axis. Due to the presence of fibers oriented than axially or
circumferentially, coupling is made between bending and twisting. The shaft has
a uniform, circular cross-section.The
following displacement field of a rotating shaft is assumed by choosing the coordinate
axisx to coincide with the shaft axis: (1)U(x,y,z,t)=U0(x,t)+zβx(x,t)−yβy(x,t),V(x,y,z,t)=V0(x,t)−zϕ(x,t),W(x,y,z,t)=W0(x,t)+yϕ(x,t), where U, V, and W are the
flexural displacements of any point on the cross-section of the shaft in the x, y, and z directions (see Figure 1), the variables U0, V0, and W0 are the flexural displacements of
the shaft's axis, while βx and βy are the rotation angles of the cross-section, about the y and z axis, respectively. The ϕ is the angular displacement of the
cross-section due to the torsion deformation of the shaft.Figure 1
The cylindrical coordinate system.The
strain components in the cylindrical coordinate system (as shown in Figures1-2) can be  written in terms of the displacement variables defined earlier as (2)εxx=∂U0∂x+rsinθ∂βx∂x−rcosθ∂βy∂x,εrr=εθθ=εrθ=0,εxθ=εθx=12(βysinθ+βxcosθ−sinθ∂V0∂x+cosθ∂W0∂x+r∂ϕ∂x),εxr=εrx=12(βxsinθ−βycosθ−sinθ∂W0∂x+cosθ∂V0∂x).Figure 2
k-layers of composite shaft.Let
us consider a composite shaft consists ofk layered (Figure 2) of fiber
inclusion reinforced laminate. The constitutive relations for each layer are
described by
(3){σxxσθθσrrτrθτxrτxθ}=[C11′C12′C13′00C16′C12′C22′C23′00C26′C13′C23′C33′00C36′000C44′C45′0000C45′C55′0C16′C26′C36′00C66′]{εxxεθθεrrγrθγxrγxθ}, where Cij′ are the
effective elastic constants, they are related to
lamination angle η (as shown in Figures 3-4) and the elastic
constants of principal axes.Figure 3
A typical composite lamina
and its principal axes.Figure 4
The definitions of the principal coordinate axes on an arbitrary layer
of the composite.The
stress-strain relations of thenth layer expressed in the cylindrical
coordinate system (Figure 5) can be expressed as
(4)σxx=C11n′εxx+ksC16n′γxθ,τxθ=τθx=ksC16n′εxx+ksC66n′γxθ,τxr=τrx=ksC55n′γxr, where ks is the
transverse shear correction factor.The stress components: (a)  in the coordinate axes (x,y,z), (b)  in the coordinate axes (x,r,θ).
(a)(b)The formula of the strain energy is(5)Ed=12∫V(σxxεxx+2τxrεxr+2τxθεxθ)dV.The various components of strain energy come from the shaft:(6)Ed=12A11∫0L(∂U0∂x)2dx+12B11[∫0L(∂βx∂x)2dx+∫0L(∂βy∂x)2dx]+12ksB66∫0L(∂ϕ∂x)2dx+12ksA16×[2∫0L∂ϕ∂x∂U0∂xdx+∫0Lβy∂βx∂xdx−∫0Lβx∂βy∂xdx−∫0L∂V0∂x∂βx∂xdx−∫0L∂W0∂x∂βy∂xdx]+12ks(A55+A66)×[∫0L(∂V0∂x)2dx+∫0L(∂W0∂x)2dx+∫0Lβx2dx+∫0Lβy2dx+2∫0Lβx∂W0∂xdx−2∫0Lβy∂V0∂xdx],
where  Aij and Bij are given in the appendix.The  kinetic energy of the rotating composite shaft, including the effects of translatory and rotary inertia, can be
written as(7)Ec=12∫0L[Im(U˙02+V˙02+W˙02)+Id(β˙x2+β˙y2)−2ΩIpβxβ˙y+2ΩIpϕ˙+Ipϕ˙2+Ω2Ip+Ω2Id(βx2+βy2)]dx, where Ω is the rotating speed of the shaft
which is assumed constant, L is the length of the shaft, the 2ΩIpβxβ˙y term accounts for the gyroscopic effect, and Id(β˙x2+β˙y2) represent the
rotary inertia effect. The mass moments of inertia Im, the
diametrical mass moments of inertia Id, and polar
mass moment of inertia Ip of rotating shaft per unit length
are defined in the appendix. As the Ω2Id(βx2+βy2) term is far smaller than Ω2Ip, it will neglected in further analysis.


## 2.2. Hierarchical Beam Element Formulation

The spinning flexible beam is descretised into one hierarchical finite element is shown
in Figure6. The element's nodal d.o.f. at each node are U0, V0, W0, βx, βy,  and ϕ. The local and nondimensional coordinates are
related by
(8)ξ=x/L,(0≤ξ≤1). The vector displacement formed by the variables U0, V0, W0, βx, βy,  and ϕ can be written as
(9)U0=∑m=1pUxm(t)⋅fm(ξ),V0=∑m=1pVym(t)⋅fm(ξ),W0=∑m=1pWzm(t)⋅fm(ξ),βx=∑m=1pβxβxm(t)⋅fm(ξ),βy=∑m=1pβyβym(t)⋅fm(ξ),ϕ=∑m=1pϕϕm(t)⋅fm(ξ). And it can be expressed as
(10){U0V0W0βxβyϕ}=[N]{q}, where
[N] is the matrix of the shape functions, given by
(11)[N]=[[NU]000000[NV]000000[NW]000000[Nβx]000000[Nβy]000000[Nϕ]],(12)[NU,V,V,βx,βy,ϕ]=[f1f2⋯fpU,pV,pW,pβx,pβy,pϕ], where pU, pV, pW, pβx, pβy,  and pϕ are the numbers of hierarchical terms of displacements (are the numbers of shape functions of displacements).Figure 6
Beam element
with two nodes.In this work,pU=pV=pW=pβx=pβy=pϕ=p.The
vector of generalised coordinates given by(13){q}={qU,qV,qW,qβx,qβy,qϕ}T, where
(14){qU}={x1,x2,x3,…,xpU}Texp(jωt),{qV}={y1,y2,y3,…,ypV}Texp(jωt),{qW}={z1,z2,z3,…,ypW}Texp(jωt),{qβx}={βx1,βx2,βx3,…,βxPβx}Texp(jωt),{qβy}={βy1,βy2,βy3,…,βyPβy}Texp(jωt),{qϕ}={ϕ1,ϕ2,ϕ3,…,ϕpϕ}Texp(jωt).The
group of the shape functions used in this study is given by(15)f1=1−ξ,f2=ξ,fr+2=sin(δrξ),δr=rπ;r=1,2,3,….The functions (f1, f2) are
those of the finite element method necessary to describe the nodal
displacements of the element, whereas the trigonometric functions fr+2 contribute only to the internal field of displacement and do not affect nodal
displacements. The most attractive particularity of the trigonometric functions
is that they offer great numerical stability. The shaft is modeled by only one
element called hierarchical finite element.By applying the Euler-Lagrange equations, the equations of
motion of free vibration of spinning flexible shaft can be obtained:(16)[M]{q¨}+[G]{q˙}+[K]{q}={0}.The system obtained is linear of which equations are coupled
once by the gyroscopic effect, represented by the matrix [G].[M] and [K] are the conventional hierarchical finite element
mass and stiffness matrix, [G] is the gyroscopic matrix (the different matrices
of the system of equation are given in the appendix).


## 3. Results

A program based on the formulation proposed
was developed for the resolution of (16).


## 3.1. Convergence

First,
the material properties for boron-epoxy are listed in Table1, and the
geometric parameters
are L=2.47m, D=12.69cm, e=1.321mm, 10 layers of
equal thickness (90°, 45°, −45°, 0∘6, 90°). The shear correction
factor ks=0.503 and the rotating speed Ω=0.Table 1
Properties
of composite materials [22].
Boron-epoxyGraphite-epoxyE11 (GPa)211.0139.0E22 (GPa)24.111.0G12 (GPa)6.96.05G23 (GPa)6.93.78v120.360.313ρ (kg/m3)1967.01578.0The results of the five bending modes for
various boundary conditions of the composite shaft as a function of the number
of hierarchical terms are shown in Figures7 and 8. Figures clearly show that
rapid convergence from above to the exact values occur as the number of hierarchical
terms is
increased.Figure 7
Convergence
of the frequencyω for the 5 bending modes of the simply-supported
(S-S) shaft as a function of the number of hierarchical terms p.Figure 8
Convergence of the frequencyω for the 5 bending
modes of the clamped-clamped (C-C) shaft as a function of the number of hierarchical
terms p.


## 3.2. Numerical Examples and Discussions

In all the following examples,p=10.In
the first example, the properties of the boron-epoxy composite shaft are given
by Table1. The results obtained using the present model are shown in Table 2
together with those of referenced papers. As can be seen from the table, our
results are close to those predicted by other beam theories. Since in the studied
example the wall of the shaft is relatively thin, models based on shell
theories [4] are expected to yield more accurate results. In the present
example, the critical speed measured from the experiment however is still
underestimated by using Sander shell theory, while overestimated by Donnell
shallow shell theory. When the material of the shaft is changed to the graphite-epoxy given in Table 1 with other conditions left unchanged, the critical speed obtained from the
present model is shown in Table 3.
In this case, the result from the present model is
compatible to that of the Bresse-Timoshenko beam theory of M.-Y. Chang et al. [6].Table 2
The critical speed of the
boron-epoxy composite shaft.
L=2.47m, D=12.69cm, e=1.321mm, 10 layers of equal thickness (90°, 45°, −45°, 0°6,
90°), ks=0.503Theory  or methodFirst  criticalspeed (rpm)Zinberg and Symonds [1]Measured experimentally6000EMBT5780dos Reis et al. [2]Bernoulli-Euler beam
theory with stiffness4942determined by shell finite elementsKim and Bert [4]Sanders shell theory5872Bert [23]Donnell shallow shell
theory6399Bert and Kim [22]Bernoulli-Euler beam
theory5919Bresse-Timoshenko beam
theory5788Gupta and Singh [3]EMBT5747LBT5620M.-Y. Chang et al. [6]Continuum-based
Timoshenko beam theory5762PresentTimoshenko beam theory by thep-version5760of the finite
element methodTable 3
The critical speed of the
graphite-epoxy composite shaft.
L=2.47m, D=12.69cm, e=1.321mm,  10 layers of equal thickness (90°, 45°, −45°, 0°6,
90°),  ks=0.503Theory
or methodFirst
criticalspeed (rpm)Bert and Kim [22]Sanders shell theory5349Donnell shallow shell
theory5805Bernoulli-Euler beam
theory5302Bresse-Timoshenko beam
theory5113M.-Y. Chang et al. [6]Continuum-based
Timoshenko beam theory5197PresentTimoshenko beam theory by thep-version5200of the finite
element methodNext,
comparisons are made with those of [6, 22] for different length to mean
diameter ratios L/D. The shafts being analysed are made of the graphite-epoxy material
given in Table 1 and all have the same lamination [90°/45°/45°/0∘6/90°].
The mean diameter and the wall thickness of the shaft remain the same as the
previous examples. The shear correction factor being used is again 0.503. The
results are listed in Table 4. Further results being compared are for
generalized orthotropic composite tube of different lamination angles  η.
The results are shown in Table 5.Table 4
The critical speed (rpm) of the
graphite-epoxy composite shaft for various lengths to mean diameter ratio.
Theory or methodL/D25101520253035Sanders shell112400416801645085855183344124401816Bert and Kim [22]Bernoulli-Euler329600768202021090725121328322821677Bresse-Timoshenko176300548301788085434945320922461658M.-Y. Chang et al. [6]Continuum-based
Timoshenko181996557061792985274925319222331648beam theoryPresentTimoshenko beam theory by184667561961800585494934319822361650thep-version of the finiteelement methodTable 5
The critical speed (rpm) of the
graphite-epoxy composite shaft for various lamination angles.
Theory or methodLamination angleη (°)0153045607590Sanders shell5527436533082386212020201997Bert and Kim [22]Bernoulli-Euler6425539342693171229218851813Bresse-Timoshenko6072520941973143227818741803M.-Y. Chang
et al. [6]Continuum-based
Timoshenko6072533142063124228418901816beam theoryTimoshenko beam theory by the6094535942223129228418901816Presentp-version of the finite
elementmethodIn
our work, the shaft is modelled by only one element with two nodes, but in the
model of [6] the shaft is modelled by 20 finite elements of equal length (h-version).
The rapid convergence while taking only one element and a reduced number of
shape functions shows the advantage of the method used.From
the thin-walled shaft systems studied above, the present shaft model yields
result in all cases close to those of the model of Bert and Kim [22] based on
the Bresse-Timoshenko theory. One should stress here that the present model is
not only applicable to the thin-walled composite shafts as studied above, but also
to the thick-walled shafts as well as to the solid ones.In
the following example, the frequencies responses of a graphite-epoxy composite
shaft system are analysed. The material properties are those listed in Table1.
The lamination scheme of the shaft remains the same as previous examples, while
its geometric properties, the Campbell
diagram containing the frequencies of the first five pairs of bending whirling
modes of the above composite system is shown in Figure 9. The intersection
point of the line (Ω = ω) with the whirling frequency curves indicate the speed
at which the shaft will vibrate violently (i.e., the critical speed). The first
10 eigenvalues correspond to 5 forward (F) and 5 backward (B) whirling bending
modes of the shaft.Figure 9
The Campbell
diagram of a composite shaft (F)
forward modes, (B) backward modes).The gyroscopic effect causes a coupling of orthogonal displacements to the axis of rotation, and by consequence separates the frequencies in two branches: backward precession mode and forward precession mode. In all cases, the backward modes increase
with increasing rotating speed however the forward modes decrease.Figure10 shows the variation of the bending fundamental frequency ω as a function of the rotating speed Ω (Campbell  diagram) for different boundary conditions and different rotating speeds.Figure 10
The first backward (1B) and forward (1F) bending mode of a composite
shaft for different boundary
conditions and different rotating speeds (S: simply supported; C: clamped; F: free).The gyroscopic effect inherent to rotating structures induces a precession motion. The backward precession modes (1B) increase with
increasing the rotating speed, however the forward precession (1F) modes decrease.


## 4. Conclusion

This paper has presented the free vibration analysis of spinning
composite shaft using thep-version,
hierarchical finite element method with trigonometric
shape functions. Results obtained using the method has been evaluated against
those available in the literature and the agreement has been found to be good.
The main conclusions have emerged from this work, these are itemised below.(1)
Monotonic and uniform convergence is found to occur as the number of
hierarchical functions is increased.(2)
The dynamic characteristics of rotating composite shaft are
influenced significantly by varying the fiber  orientation, the rotating speed,
and boundary conditions.(3)
The gyroscopic effect causes a coupling of orthogonal displacements to the axis of rotation, and by consequence separates the frequencies in two branches: backward and forward precession modes. In all cases,
the backward modes increase with increasing rotating speed however the forward
modes decrease.



---

*Source: 752062-2008-06-22.xml*