# Stability analysis of rotating composite shafts considering internal damping and coupling effects

Thouraya Baranger, Eric Chatelet, Marie-Ange Andrianoely, Georges

Jacquet-Richardet

► To cite this version:

Thouraya Baranger, Eric Chatelet, Marie-Ange Andrianoely, Georges Jacquet-Richardet. Stability analysis of rotating composite shafts considering internal damping and coupling effects. 8th IFTOMM International Conference on Rotordynamics, Sep 2010, seoul, South Korea. pp.1-8. hal-04778919

HAL Id: hal-04778919

https://hal.science/hal-04778919v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# STABILITY ANALYSIS OF ROTATING COMPOSITE SHAFTS CONSIDERING INTERNAL DAMPING AND COUPLING EFFECTS

**Thouraya Nouri Baranger**

Université de Lyon, CNRS
LaMCoS UMR5259, Université Lyon 1
F-69622 Villeurbanne, France

**Marie-Ange Andrianoely**

Université de Lyon, CNRS
LaMCoS UMR5259, INSA-Lyon,
F-69621 Villeurbanne, France

**Eric Chatelet**

Université de Lyon, CNRS
LaMCoS UMR5259, INSA-Lyon,
F-69621 Villeurbanne, France

**Georges Jacquet-Richardet**

Université de Lyon, CNRS
LaMCoS UMR5259, INSA-Lyon,
F-69621 Villeurbanne, France

## ABSTRACT

This paper is concerned with the dynamic instability of an internally damped rotating composite shaft. A tridimensional homogenized finite element beam model, taking into account both internal damping and coupling effects, is developed and, then used to evaluate natural frequencies and instability thresholds of the rotating structure. The influence of laminate parameters: stacking sequences, fiber orientation on the natural frequencies and the instability thresholds of the shaft is analyzed. First, these parametric studies show that shaft instability thresholds can be very sensitive to laminate parameters. Then the influence of coupling effects on the instability thresholds is investigated.

## KEYWORDS

Rotordynamics, Composite, Internal Damping, Stability, Shear and Coupling effects

## INTRODUCTION

Composite materials have interesting properties such as high strength to weight ratio, compared to metals, which make them very attractive for rotating systems. Attempts are being made to replace metal shafts by composite ones in many applications: drive shafts for helicopters, centrifugal separators, and cylindrical tubes for the automotive and marine industries (Zorzi and Giordano [1], Darlow [2], Gupta and Singh [3, 4, 5], Chatelet [6]). They also provide designers with the possibility of obtaining predetermined behaviours, in terms of position of critical speeds, by changing the arrangement of the different composite layers: orientation and number of plies (Bauchau [7], Gubran and Gupta [8], Chatelet [6] and Pereira [9]).

On the other hand, these materials have relatively high-damping characteristics. For a rotor made with composite materials, internal damping is much more significant than when associated with a metal rotor. Unfortunately, such damping may cause instability as shown by Wettergren [10].

Accurate prediction of damping characteristics of rotor systems is therefore fundamental in the design of rotating machines as it provides estimations on safe-ranges of speeds of rotation. Over the last few years, many studies have focused on predicting critical speeds, natural frequencies, unbalance responses and, in particular, instability thresholds. Newkirk [11] observed that rotor-disk systems would undergo violent whirling at the first natural frequency at speeds above the first critical speed. Kimball [12] showed that internal damping destabilizes the whirling motion if the rotation speed of the rotor exceeds the first critical speed. In addition, Bucciarelli [13] showed that the instability criterion based on the ratio of energy dissipated between internal and external damping is inaccurate and that internal forces can produce instability by coupling spin and whirl motions.

Classical results have been obtained and showed that rotor stability is improved by increasing external damping, whereas increasing internal damping may reduce the instability threshold. However, most of the published studies deal with metal rotating structures and remain exclusively numerical without precise estimations of internal damping.

Several finite element formulations have been performed for the analysis of composite shafts. These formulations are based on homogenized beam and shell theories. The equivalent modulus beam theory (EMBT), which is widely used for the dynamic analysis of composite shafts, was firstly introduced by Tsai [14]. With this approach, equivalent longitudinal Young and in-plane shear moduli are identified by using laminate theory for symmetrical stacking. Then, classical beam theory can be used to model the shaft, see Pereira [9] and Singh and Gupta [4]. This approach has many limitations which are summarized by Singh and Gupta in [3]. They studied the natural frequencies and damping ratios in flexural modes of cylindrical laminate tubes and compared shell and EMBT models for symmetric laminate stacking, concluding that in the case of the tube configurations usually used in composite shaft applications, the differences in flexural frequencies between the two models are negligibly small.
Using shell theory in [15], the same authors, showed that the modal loss factors are more sensitive to parametric (laminate stacking, angle orientations, etc.) changes than frequency values. They also presented in [5] a comparison between EMBT theory and Layerwise Beam Theory (LBT) for symmetric and asymmetric stacking. They showed that LBT is more efficient that EMBT because it takes into account the effect of changed stacking, thickness shear deformation and bending-stretching coupling. However, LBT requires the development of a complex beam element with a high number of degrees of freedom dependent on the number of layers, making the method too expensive. Recently, Gubran and Gupta presented in [8] a modified EMBT method which takes into account the effects of a stacking sequence and different coupling mechanisms. They considered a Graphite/Epoxy shaft simply supported on rigid bearings and compared the first three frequencies with those obtained by using the LBT method. In spite of its simplicity, the natural frequencies obtained using modified EMBT excluding different coupling effects agree well with those obtained using LBT and those reported in the literature. In these cited works, the internal damping is not often taken into account, except in [3] where viscoelastic material damping is assumed.

In this paper, a three dimensional homogenized beam is introduced in order to be used to analysis frequencies and instability thresholds. In this approach, only beam theory assumptions are considered. Then coupling effects induced by the fiber orientations appears in the equation of equilibrium. These coupling effects concern: longitudinal and torsion deformation and out of plane bending. In this approach, the orthotropic effect of the internal damping of each ply is taken into account.

In the following, an outline of the formulation is presented. In section 2, equations of motion of rotordynamics with and without internal damping are presented and compared. In section 3, the orthotropic properties of a layer component of a composite shaft are presented, and are expressed in the beam coordinate system. A homogenized beam theory is developed; elastic energy and dissipative virtual work are given. Then, in section 4, the homogenized beam parameters are expressed as a function of the layer parameters. Numerical applications are presented in section 5. In the first one, the approach is validated via an isotropic metallic shaft presented in Lalanne [16]. The second application is that presented in the work of Pereira [9]. Shaft frequencies and instability thresholds regarding stacking order and fiber orientations are outlined. These results are also compared to those obtained from a classical equivalent modulus beam theory. The third application considering a symmetrical and balanced configuration illustrates the influence of coupling effects on the estimation of the instability threshold.

## ROTORDYNAMICS

In the fixed frame, the following equations of motion are associated with a rotor made of an isotropic material (Lalanne [16]):

$$[M]\{\ddot{d}\} + [C + G(\Omega)]\{\dot{d}\} + [K]\{d\} = \{F(t)\} \quad (1)$$

Where [M] is the symmetric mass matrix, [C] the external damping matrix, [G(Ω)] is the global asymmetric matrix including an anti-symmetric gyroscopic matrix (function of Ω speed of rotation) and a frequently asymmetric matrix due to the characteristics of the bearings, [K] is the elastic stiffness matrix that is frequently asymmetric due to the characteristics of bearings, {F(t)} is the generalized force vector {dd}, {d} and {d} are respectively nodal acceleration, velocity and displacement vectors. Taking into account the material's dissipative effects gives two other matrices associated with internal damping, as shown in Sino [17]:

$$[M]\{\ddot{d}\} + [C_i + C + G(\Omega)]\{\dot{d}\} + [K + K_i(\Omega)]\{d\} = \{F(t)\} \quad (2)$$

where [C_i] is the internal damping matrix and [K_i(Ω)] a stiffness matrix which depends on the internal damping and also on the rotational speed Ω. The anisotropic properties of composite materials and their lightness can be used to optimize composite shafts in order to improve their dynamic behavior. Compared to metals, composite materials have higher damping capacities which can induce a destabilizing effect on the rotor motion. When modeling composite rotors with the equivalent modulus beam theory, equation is considered directly. However, as mentioned earlier, this approach is based on symmetric stacking laminate theory and cannot take into account the influence of layer stacking order. In the following, a more general homogenized beam model is proposed.

## COMPOSITE ROTOR

The shaft studied can be obtained by wounding several plies of embedded fibbers on a mandrel, as shown in figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Composite Rotor</figcaption></figure>

Each ply has an orthotropic mechanical behavior, as shown in figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: Plan of ply</figcaption></figure>
The generalized Hooke's law for an orthotropic material is written as follows:

$$ \{\sigma\} = [\mathcal{Q}]\{\epsilon\} \quad \text{or} \quad \{\epsilon\} = [\mathcal{S}]\{\sigma\} \qquad (3) $$

where $\{\sigma\}$ and $\{\epsilon\}$ are respectively the stress and strain fields, $[\mathcal{Q}]$ and $[\mathcal{S}]$ are respectively the material stiffness and compliance matrices. Only the expression of the compliance matrix will be developed here. The stiffness matrix can be obtained by considering that $[\mathcal{Q}] = [\mathcal{S}]^{-1}$. When linked to the orthotropic axis, Hooke's law takes the following form:

$$ \begin{bmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \gamma^{23} \\ \gamma^{13} \\ \gamma^{12} \end{bmatrix} = \begin{bmatrix} \frac{1}{E_1} & -\frac{\nu_{21}}{E_2} & -\frac{\nu_{31}}{E_3} & 0 & 0 & 0 \\ -\frac{\nu_{12}}{E_1} & \frac{1}{E_2} & -\frac{\nu_{32}}{E_3} & 0 & 0 & 0 \\ -\frac{\nu_{13}}{E_1} & -\frac{\nu_{23}}{E_2} & \frac{1}{E_3} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{23}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{13}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{12}} \end{bmatrix} \begin{bmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma^{23} \\ \sigma^{13} \\ \sigma^{12} \end{bmatrix} \quad (4) $$

where (1, 2, 3) are the orthotropic axes. 1 is the fiber direction, 2 is the direction transversal to the fibers in the ply and 3 is the perpendicular direction to the ply. The fiber direction 1 makes an angle $\Phi$ with respect to the x-axis. Considering transversally isotropic material, the following parameters have to be identified for each ply: $E_1$ and $E_2 = E_3$ Young moduli in the orthotropic axes; $G_{23}$, $G_{13} = G_{12}$, transversal shear moduli, and $\nu_{21} = \nu_{13}$ Poisson's ratios.

The behavior of a viscoelastic composite material in harmonic steady-state motion can be described by the complex constitutive relation. Assuming cyclic loading, the complex stress component is written as:

$$ \{\sigma\} = [\mathcal{Q}]\{\epsilon\} + j[\mathcal{Q}]\{\dot{\epsilon}\} \qquad (5) $$

with $[\mathcal{Q}] = [\eta][\mathcal{Q}]$, $[\eta]$ is the damping matrix of the ply and $j$ is the imaginary unit. The dissipative properties of a ply can also be expressed by using the specific damping capacities matrix $[\Psi]$. Usually energy dissipation in solids is characterized by the relative energy dissipation which is defined as the ratio of the energy losses $\Delta W$ in a unit volume of a body, to the elastic energy $W$ under a given stress-strain state, Zinoviev [18]:

$$ \psi = \frac{\Delta W}{W} = \frac{\int_{0}^{2\pi/\alpha} \{\epsilon\}^{t} [S^{\prime}][\epsilon] dt}{\int_{0}^{2\pi/\alpha} \{\epsilon\}^{t} [S^{\prime}][\epsilon] dt} \qquad (6) $$

where $[S^{\prime}]$ is the damped compliance matrix. The composite ply has three main directions of specific damping capacity, which can be expressed by the following matrix.

$$ [\Psi] = \begin{bmatrix} \psi_{11} & 0 & 0 & 0 & 0 & 0 \\ 0 & \psi_{22} & 0 & 0 & 0 & 0 \\ 0 & 0 & \psi_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & \psi_{23} & 0 & 0 \\ 0 & 0 & 0 & 0 & \psi_{13} & 0 \\ 0 & 0 & 0 & 0 & 0 & \psi_{12} \end{bmatrix} \qquad (7) $$

Here $\psi_{ii}$ are the specific damping capacities associated longitudinal effects and $\psi_{ij}$ are those associated with the transversal shear effect. These coefficients can be identified experimentally. The damping matrix $[\eta]$ can be linked to the specific damping capacity as follows:

$$ [\eta] = \frac{1}{2\pi} [\Psi] \qquad (8) $$

Consequently, the damped material stiffness matrix $[\mathcal{Q}^{\prime}]$ is expressed as a function of the specific damping capacity as follows:

$$ [\mathcal{Q}^{\prime}] = \frac{1}{2\pi} [\Psi] [\mathcal{Q}] \qquad (9) $$

Consider the cylindrical coordinate system $(x, r, \theta)$ as mentioned figure 1. Each ply $p$ is characterized by the angle $\Phi_p$ between the shaft axis $x$ and the fiber axis $l$ axes. The stress-strain relation in this chosen cylindrical coordinate system can be written as:

$$ [\bar{Q}] = [T][\mathcal{Q}][T]^T \qquad (10) $$

Where the transfer matrix is given by:

$$ [T] = \begin{bmatrix} c^2 & s^2 & 0 & 0 & 0 & 2cs \\ s^2 & c^2 & 0 & 0 & 0 & -2cs \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & c & s & 0 \\ 0 & 0 & 0 & -s & c & 0 \\ -cs & cs & 0 & 0 & 0 & (c^2 - s^2) \end{bmatrix} \qquad (11) $$

with $c = \cos(\Phi_p)$ and $s = \sin(\Phi_p)$. Then, the stiffness matrix takes the following form and coupling terms appear.

$$ \begin{Bmatrix} \sigma_{xx} \\ \sigma_{rr} \\ \sigma_{\theta\theta} \\ \sigma_{r\theta} \\ \sigma_{xr} \\ \sigma_{x\theta} \end{Bmatrix} = \begin{bmatrix} \bar{Q}_{11} & \bar{Q}_{12} & \bar{Q}_{13} & 0 & 0 & \bar{Q}_{16} \\ \bar{Q}_{12} & \bar{Q}_{22} & \bar{Q}_{23} & 0 & 0 & \bar{Q}_{26} \\ \bar{Q}_{13} & \bar{Q}_{23} & \bar{Q}_{33} & 0 & 0 & \bar{Q}_{36} \\ 0 & 0 & 0 & \bar{Q}_{44} & \bar{Q}_{45} & 0 \\ 0 & 0 & 0 & \bar{Q}_{45} & \bar{Q}_{55} & 0 \\ \bar{Q}_{16} & \bar{Q}_{26} & \bar{Q}_{36} & 0 & 0 & \bar{Q}_{66} \end{bmatrix} \begin{Bmatrix} \epsilon_{xx} \\ \epsilon_{rr} \\ \epsilon_{\theta\theta} \\ \epsilon_{r\theta} \\ \epsilon_{xr} \\ \epsilon_{x\theta} \end{Bmatrix} \qquad (12) $$

Let's consider a multilayered composite shaft made of $N$ orthotropic layers. If the stacking sequence is symmetric the shaft has a typical beam behavior and can be modeled by using Timoshenko beam theory. If the stacking sequence is nonsymmetric, mechanical coupling effects such as bending-stretching, twisting-stretching and shear-stretching will be present. In this cylindrical coordinate system, the beam theory assumption leads to the following equality:

$$ \sigma_{rr} = \sigma_{r\theta\theta} = \sigma_{\theta\theta} = 0 \qquad (13) $$

Thus, the stress strain relation can be written as follows:

$$ \begin{Bmatrix} \sigma_{xx} \\ \sigma_{x\theta} \\ \sigma_{xr} \end{Bmatrix} = \begin{bmatrix} \tilde{Q}_{11} & \tilde{Q}_{16} & 0 \\ \tilde{Q}_{16} & \tilde{Q}_{66} & 0 \\ 0 & 0 & \tilde{Q}_{55} \end{bmatrix} \begin{Bmatrix} \epsilon_{xx} \\ \epsilon_{x\theta} \\ \epsilon_{xr} \end{Bmatrix} \qquad (14) $$

The components $\tilde{Q}_{ij}$ are easily obtained from the equation (11) by considering the above assumption. The
relation (12) can be written in the Cartesian coordinate system of the shaft as follows:

$$
\left\{ \begin{matrix} \sigma_{xx} \\ \sigma_{xy} \\ \sigma_{xz} \end{matrix} \right\} = \frac{\begin{bmatrix} \tilde{Q}_{11} & n\tilde{Q}_{16} & -n\tilde{Q}_{16} \\ n\tilde{Q}_{16} & m^2\tilde{Q}_{55} + n^2\tilde{Q}_{66} & mn(-\tilde{Q}_{66} + \tilde{Q}_{55}) \\ -m\tilde{Q}_{16} & mn(-\tilde{Q}_{66} + \tilde{Q}_{55}) & m^2\tilde{Q}_{55} + n^2\tilde{Q}_{66} \end{bmatrix}}{[c]} \left\{ \begin{matrix} \epsilon_{xx} \\ \epsilon_{xy} \\ \epsilon_{xz} \end{matrix} \right\} \quad (15)
$$

With $m=\cos(\theta)$, $m=\sin(\theta)$, and the angle $\theta$ is the second coordinate of the cylindrical coordinate system.

Let's consider a beam theory to model the composite rotor illustrated in figure 3. Thus the shaft is modeled as a beam with a constant circular cross-section. The finite element considered here has 6 degrees of freedom for each node: three translations $u, v$ and $w$ and three rotations $\theta_k, \theta_l$ and $\theta_m$. The continuous displacement field at material points along the rotor cross-section is described as follows:

$$
\{u(x, y, z)\} = \begin{cases} u_x = u + z\theta_y - y\theta_z \\ u_y = v - z\theta_x \\ u_z = w + y\theta_x \end{cases} \quad (16)
$$

Hence, the deformation field has the following form:

$$
\{\epsilon\} = \begin{cases} \epsilon_{xx} = \frac{\partial u}{\partial x} - y \frac{\partial \theta_z}{\partial x} + z \frac{\partial \theta_y}{\partial x} \\ \gamma_{xz} = \theta_y + \frac{\partial w}{\partial x} + y \frac{\partial \theta_x}{\partial x} \\ \gamma_{xy} = -\theta_z + \frac{\partial v}{\partial x} - z \frac{\partial \theta_x}{\partial x} \end{cases} \quad (17)
$$

Consequently, for each ply $k$ of the rotor cross-section, the stress-strain relation is written as:

$$
\{\sigma\} = [C]\{\epsilon\} + [C"]\{\dot{\epsilon}\} \quad (18)
$$

$$
\text{Where: } [C''] = [\eta][C] \qquad (19)
$$

The virtual elastic and dissipative work has the following expression:

$$
\delta W = \int_0^L \int_S (\sigma_{xx} \delta \epsilon_{xx} + \sigma_{xz} \delta \gamma_{xz} + \sigma_{xy} \delta \gamma_{xy}) dSdx \quad (20)
$$

where S is the cross section.

HOMOGENIZATION

The rotor has constant geometric properties along its
longitudinal x-axis. The homogenized mechanical
characteristics are extracted from equations (16) to (20) by
evaluating the integrals over the cross-section. We denote by
{ε⁰} the deformations associated to the x-axis such that:

$$
\{\epsilon^0\}^T = \left\langle \frac{\partial u}{\partial x} \quad \frac{\partial v}{\partial x} - \theta_z \quad \frac{\partial w}{\partial x} + \theta_y \quad \frac{\partial \theta_x}{\partial x} \quad \frac{\partial \theta_y}{\partial x} \quad \frac{\partial \theta_z}{\partial x} \right\rangle \quad (21)
$$

And by {$F^0$} the resulting forces and moments acting on the
x-axis such that:

$$
\{F^0\}^T = \langle N \ T_y \ T_z \ M_x \ M_y \ M_z \rangle \quad (22)
$$

Then the elastic and dissipative energies can be written as
follows:

$$
\delta W = \langle F^0 \rangle ([C_e] + \Omega[C_{a1}])\{\epsilon^0\} + \langle F^0 \rangle ([C_{a2}])\{\dot{\epsilon}^0\} \quad (23)
$$

Where:

$$
[\tilde{C}_e] = \begin{bmatrix} \tilde{C}_{11}^e & 0 & 0 & \tilde{C}_{14}^e & 0 & 0 \\ 0 & \tilde{C}_{22}^e & 0 & 0 & \tilde{C}_{25}^e & 0 \\ 0 & 0 & \tilde{C}_{22}^e & 0 & 0 & \tilde{C}_{25}^e \\ \tilde{C}_{14}^e & 0 & 0 & \tilde{C}_{44}^e & 0 & 0 \\ 0 & \tilde{C}_{25}^e & 0 & 0 & \tilde{C}_{55}^e & 0 \\ 0 & 0 & \tilde{C}_{25}^e & 0 & 0 & \tilde{C}_{55}^e \end{bmatrix} \quad (24)
$$

$$
[\tilde{C}_{a1}] = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \tilde{C}_{22}^{a1} & 0 & 0 & \tilde{C}_{26}^{a1} \\ 0 & -\tilde{C}_{22}^{a1} & 0 & 0 & -\tilde{C}_{26}^{a1} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \tilde{C}_{26}^{a1} & 0 & 0 & \tilde{C}_{56}^{a1} \\ 0 & -\tilde{C}_{26}^{a1} & 0 & 0 & -\tilde{C}_{56}^{a1} & 0 \end{bmatrix} \quad (25)
$$

$$
[\tilde{C}_{a2}] = \begin{bmatrix} \tilde{C}_{11}^{a2} & 0 & 0 & \tilde{C}_{14}^{a2} & 0 & 0 \\ 0 & \tilde{C}_{22}^{a2} & 0 & 0 & \tilde{C}_{25}^{a2} & 0 \\ 0 & 0 & \tilde{C}_{22}^{a2} & 0 & 0 & \tilde{C}_{25}^{a2} \\ \tilde{C}_{14}^{a2} & 0 & 0 & \tilde{C}_{44}^{a2} & 0 & 0 \\ 0 & \tilde{C}_{25}^{a2} & 0 & 0 & \tilde{C}_{55}^{a2} & 0 \\ 0 & 0 & \tilde{C}_{25}^{a2} & 0 & 0 & \tilde{C}_{55}^{a2} \end{bmatrix} \quad (26)
$$

With:

$$
\begin{align}
\tilde{C}_{11}^e &= \pi \sum_{k=1}^{p} C_{11}^k (r_{2k}^2 - r_{1k}^2) \notag \\
\tilde{C}_{14}^e &= -\frac{2\pi}{3} \sum_{k=1}^{p} C_{12}^k (r_{2k}^3 - r_{1k}^3) \notag \\
\tilde{C}_{25}^e &= \frac{\pi}{3} \sum_{k=1}^{p} C_{12}^k (r_{2k}^3 - r_{1k}^3) \tag{27}
\end{align}
$$

$$
\tilde{C}_{22}^{e} = \frac{\pi}{2} \sum_{k=1}^{p} (C_{22}^{k} + C_{33}^{k})(r_{2k}^{2} - r_{1k}^{2})
$$

$$
\tilde{C}_{44}^e = \pi \sum_{k=1}^{p} C_{22}^k (r_{2k}^4 - r_{1k}^4)
$$

$$
\tilde{C}_{55}^{e} = \frac{\pi}{4} \sum_{k=1}^{p} C_{11}^{k} (r_{2k}^{4} - r_{1k}^{4})
$$

Where $C_{ij}^k$ are the coefficient of the matrix $[C]$ for the $k^{th}$ ply.
For $[\tilde{C}_{a2}]$, the terms have the same expression expect that
the coefficients $C_{ij}^k$ are replaced by those of $[C'']$. The terms
of $[\tilde{C}_{a1}]$ are:

$$
\tilde{C}_{23}^{a1} = \frac{\pi}{2} \sum_{k=1}^{p} (C_{22}^{*k} + C_{33}^{*k})(r_{2k}^2 - r_{1k}^2) \quad (28)
$$

$$
\tilde{C}_{26}^{a1} = \frac{\pi}{3} \sum_{k=1}^{p} C_{12}^{*k} (r_{2k}^3 - r_{1k}^3)
$$

$$
\tilde{C}_{56}^{a1} = \frac{\pi}{4} \sum_{k=1}^{p} C_{11}^{*k} (r_{2k}^4 - r_{1k}^4)
$$
The matrix $[\tilde{C}_e]$ expresses the constitutive homogenized stiffness of the shaft. Notice the existences of coupling coefficients $C_{ij}^{ek}$, which vanish when the materials are isotropic or when the stacking sequences of the shaft is symmetric. The matrix $[C_{a1}]$ expresses the damped stiffness; it is antisymmetric and will induce instability in the shaft. The matrix $[C_{a2}]$ expresses the damping effect.

Using these matrices in the finite element formulation we derive the matrices $[K]$, $[K_i]$ and $[C_i]$ in the formulas (2). In this paper we cannot details the finite element formulation, it will be included in a forthcoming paper.

# APPLICATIONS

## Metallic Shaft with 3 discs

The objective of this paragraph is to validate the proposed method with a shaft made of identical metallic plies. The rotor, shown Figure 4, is described in Lalanne [16]).

<figure><img src="image_5.png" /><figcaption>Figure 4: Rotor model with three rigid disks</figcaption></figure>

The rotor uses 13 shaft finite elements of the same length. The data are:

$$L_1 = 0.2m, L_2 = 0.3m, L_3 = 0.5m, L_4 = 0.3m.$$

The shaft cross-sectional radius is 0.05m. The discs are given in Table 1.

<table><caption>Table 1: Discs Data</caption><thead><tr><th>Disks</th><th>D<sub>1</sub></th><th>D<sub>2</sub></th><th>D<sub>3</sub></th></tr></thead><tbody><tr><td>Thickness m</td><td>0.05</td><td>0.05</td><td>0.06</td></tr><tr><td>Inner radius m</td><td>0.05</td><td>0.05</td><td>0.05</td></tr><tr><td>Outer radius m</td><td>0.12</td><td>0.2</td><td>0.2</td></tr></tbody></table>

The discs and the shaft are made of steel (E = 2e11 N/m², ρ = 7800 kg/m³, ν = 0.3). To respect isotropy of steel, the rotor is constituted by 10 identical plies of same thickness, with fiber orientation $Φ_0$.

The two bearings are assumed to be similar and characterized by:

$$k_{yy}=5 10^7 \text{ N/m}, k_{zz}=7 10^7 \text{ N/m},$$

$$c_{yy} = 5 10^2 \text{ N/m/s}, c_{zz} = 7 10^2 \text{ N/m/s}.$$

The Figure 5 presents the Campbell diagram obtained thanks to the three dimensional proposed method. This graph shows the evolution of the thirteen first frequencies as function of the speed of rotation. As we can see, gyroscopic effects are only noticeable for the bending modes (solid

lines). Frequencies of Torsion (-. Lines) or Longitudinal (... lines) modes stay constant whatever the rotation speed Ω. Only the first order line corresponding to the first engine order is represented.

<figure><img src="image_6.png" /><figcaption>Figure 5: Campbell Diagram</figcaption></figure>

The Table 2 compares the frequencies calculated with a one dimensional in-house code Rotorinsa® [19] with those obtained with the 3D method at 25 000 rpm. Here FW means Forward Whirl; BW backward Whirl. Contrary to Rotorinsa® (considering only a beam element with only four degrees of freedom (2 translations and two rotations)), the proposed method allows to identify torsion and traction/compression modes.

<table><caption>Table 2: Identified Frequencies at 25000 rpm.</caption><thead><tr><th>frequencies</th><th>Nature</th><th>Rotorinsa</th><th>3D Model</th><th>ε(%)</th></tr></thead><tbody><tr><td>1</td><td>1BW</td><td>55.41</td><td>55.67</td><td>0.47</td></tr><tr><td>2</td><td>1FW</td><td>67.19</td><td>67.01</td><td>0.27</td></tr><tr><td>3</td><td>2BW</td><td>157.90</td><td>158.45</td><td>0.34</td></tr><tr><td>4</td><td>2FW</td><td>193.64</td><td>193.45</td><td>0.10</td></tr><tr><td>5</td><td>3BW</td><td>249.87</td><td>251.07</td><td>0.47</td></tr><tr><td>6</td><td>1T</td><td>-</td><td>254.01</td><td>-</td></tr><tr><td>7</td><td>3FW</td><td>407.53</td><td>406.31</td><td>0.29</td></tr><tr><td>8</td><td>4BW</td><td>446.85</td><td>450.53</td><td>0.83</td></tr><tr><td>9</td><td>5BW</td><td>623.11</td><td>627.42</td><td>-</td></tr><tr><td>10</td><td>1L</td><td>-</td><td>636.73</td><td>0.69</td></tr><tr><td>11</td><td>4FW</td><td>715.22</td><td>710.14</td><td>0.71</td></tr><tr><td>12</td><td>2T</td><td>-</td><td>712.84</td><td>-</td></tr></tbody></table>

The example illustrates with very good accuracy the capabilities of the method. Now, it can be applied to a real composite structure.

## Composite Shaft with 2 rigid discs

In rotordynamic analysis, in order to emphasize the influence of internal damping and coupling effects on the instability thresholds, the Campbell diagram and instability thresholds are determined for a rotor made of a wounding shaft for any configuration of the stacking sequence. The structure, proposed by Pereira [9], is a composite shaft with two rigid steel rigid disks supported by two bearings at the ends as represented in Figure 6.
* Rotor: L=1.2m, D₀=0.096m, e=0.008m

* Disc: Dᵢ=0.096m, D₀=0.3m, h=0.05m

* Composite: 8 layers carbon/epoxy

<figure><img src="image_7.png" /><figcaption>Figure 6: Rotor in wounding shaft with two rigid disks</figcaption></figure>

Anisotropic bearings stiffness characteristics are described in Table 3 (considered without external damping). The material proprieties of each ply (carbon/epoxy) are summarized in Table 4

<table><caption>Table 3: Stiffness data of the anisotropic bearings</caption><thead><tr><th></th><th>K<sub>xx</sub><br>(N/m)</th><th>K<sub>zz</sub><br>(N/m)</th><th>K<sub>xz</sub><br>(N/m)</th><th>K<sub>zx</sub><br>(N/m)</th></tr></thead><tbody><tr><th>Anisotropic<br>bearings</th><td>1.10<sup>7</sup></td><td>1.10<sup>8</sup></td><td>0</td><td>0</td></tr></tbody></table>

<table><caption>Table 4: Material data of the shaft</caption><thead><tr><th>Material</th><th>E<sub>1</sub><br>(GPa)</th><th>E<sub>2</sub><br>(GPa)</th><th>G<sub>12</sub><br>(GPa)</th><th>ν<sub>12</sub></th><th>ρ<br>(kg/m)</th><th>ψ<sub>1</sub><br>(%)</th><th>ψ<sub>1</sub><br>(%)</th><th>ψ<sub>12</sub><br>(%)</th></tr></thead><tbody><tr><th>Carbon/<br>Epoxy</th><td>172.7</td><td>7.20</td><td>3.76</td><td>0.3</td><td>1446.2</td><td>0.45</td><td>4.22</td><td>7.05</td></tr></tbody></table>

### Validation of the method without couplings effects configuration [± 75°],

The first studied composite rotor configuration is [± 75°]₈. The results obtained by the proposed method are compared with those calculated by an equivalent modulus beam theory (EBMT) developed by Tsai [14]. Here no couplings effects are taken into account as both two other methods don't consider them. Table 5 compares the two first frequencies of bending calculated by the 3 methods.

<table><caption>Table 5: Mechanical characteristics of the shaft and Results</caption><thead><tr><th>Stacking<br>sequence<br>[± 75°]₈.</th><th>1BW<br>(Hz)</th><th>1FW<br>(Hz)</th><th>2BW<br>(Hz)</th><th>2FW<br>(Hz)</th><th>Instability<br>threshold<br>(rpm)</th></tr></thead><tbody><tr><th>EBMT</th><td>17.2</td><td>17.4</td><td>66.3</td><td>72.3</td><td>1020</td></tr><tr><th>SINO [17]</th><td>16.4</td><td>16.6</td><td>59.4</td><td>60.7</td><td>1110</td></tr><tr><th>3D method</th><td>16.1</td><td>16.6</td><td>60.95</td><td>62.3</td><td>1003</td></tr></tbody></table>

Such results are in perfect agreement with those obtained by Sino[17] and Pereira [9]. The Campbell diagram shown in Figure 7 illustrates instability that occurs just after the second critical speed.

<figure><img src="image_8.png" /><figcaption>Figure 7: Campbell diagram and instability regions for a laminate θ = 75° with anisotropic bearings</figcaption></figure>

### Validation of the method with couplings effects

The second composite rotor, in a balanced and symmetrical configuration [90°, 45, 0]s was studied thanks to the one dimensional homogenized finite element beam model developed by Sino[17]. As observed, the coupling effects have not a major influence on the first bending modes even if they generate apparition of other types of modes (longitudinal and torsion modes).

<figure><img src="image_9.png" /><figcaption>Figure 8: Campbell diagram and instability regions (noted by +) with anisotropic bearings [90°,45,0]s</figcaption></figure>

Figure 8 presents the associated Diagram Campbell; torsion modes are symbolized by “-.” lines. Bending modes are symbolized by “- -” lines and instability is associated with “+” symbol.

The table 6 lists the bending frequencies and instability thresholds calculated without and with internal damping with coupling effects taken into account. For information, instability thresholds is determined thanks the real part of complex eigenvalues. For comparison, results calculated by Sino are given too.
<table><caption>**Table 6: Frequencies and instability threshold for a balanced and symmetrical configuration [90<sub>2</sub>,45,0]<sub>s</sub>**</caption><thead><tr><th>Stacking sequence<br>[90<sub>2</sub>,45,0]<sub>s</sub></th><th>1BW<br>(Hz)</th><th>1FW<br>(Hz)</th><th>Instability<br>threshold<br>(rpm)</th></tr></thead><tbody><tr><td>SINO [17]<br>with internal damping<br>without coupling effects</td><td>37.28</td><td>39.94</td><td>5750</td></tr><tr><td>3D method<br>without internal damping<br>with coupling effects</td><td>38.78</td><td>41.81</td><td>stable</td></tr><tr><td>3D method<br>with internal damping<br>with coupling effects</td><td>40.21</td><td>43.14</td><td>2900</td></tr></tbody></table>

As we can see, the coupling effects have a major influence on the instability threshold as its value diminishes about 50 %. Finally, Figure 9 shows the shape mode associated to the first bending mode. It can be noticeable the influence of coupling effects as the mode is not anymore in a single plane, but have components out of plane.

1st Bending mode shape - Influence of coupling effects

<figure><img src="image_10.png" /><figcaption>**Figure 9: Out of plane bending mode due to coupling effects in configuration [90<sub>2</sub>,45,0]<sub>s</sub>**</figcaption></figure>

## CONCLUSIONS

This work deals with the stability analysis of an internally damped rotating composite shaft. A three dimensional homogenized beam is developed and compared to the classical Equivalent Beam Modulus Theory (EMBT) and the Simplified Homogenized Beam Theory developed in [17]. The method developed avoids the main drawbacks associated with EMBT formulation that considers only symmetrical and balanced stacking sequences and does not take into account the distance of composite layers from the neutral axis. It also takes into account internal damping by using the specific damping capacity of each ply of the composite assembly, and coupling effects induced by the fiber orientation.

First this approach is validated via an isotropic metallic shaft. A second application concerns a composite Shaft. Two different stacking sequences are analysed. In the first one, the critical speeds as well as the instability threshold obtained by the 3D method developed here are in good

agreement with those obtained in literature (Equivalent mudolus beam theory and simplified homogenized theory).

frequencies and instability thresholds regarding stacking order and fiber orientations are outlined. These results are also compared to those obtained from a classical equivalent modulus beam theory. The third application considering a symmetrical and balanced configuration illustrates the influence of coupling effects on the estimation of the instability threshold.

The study highlights that EMBT simplifications may lead to significant discrepancies in terms of frequencies. These discrepancies appear to be greater for instability thresholds. A qualitative study of the effects of various parameters on frequencies and instability thresholds was carried out. The analysis shows that although transversal shear has a minor influence on the first frequencies, its effect is much more significant for the following ones, thereby directly influencing instability thresholds.

However, this method requires some improvements to take account of the coupling effects induced by nonsymmetrical stacking. These improvements will be the subject of a forthcoming paper.

## REFERENCES

1. Zorzi, E.S. and Giordano, J.C., 1985, Composite Shaft Rotordynamic Evaluation, *The American Society of Mechanical Engineers*, Vol. 85-det-114.

2. Darlow M.S. and Creonte J., 1995, Optimal design of composite helicopter power transmission shafts with axially varying fiber lay-up, Journal of the American Helicopters society, vol. 40(2), pp. 50-56.

3. Singh S.P. and Gupta, K., 1994, Free damped flexural vibration analysis of composite cylindrical tubes using beam and shell theories, Journal of Sound and Vibration, Vol. 172(2), pp. 171-190.

4. Singh S.P. and Gupta K., 1996, Dynamic Analysis of Composite rotors, International Journal of Rotating Machinery, vol. 2(3), pp. 179-186.

5. Singh S.P. and Gupta K., 1996, Composite shaft rotordynamic analysis using a layerwise theory, Journal of Sound and Vibration, vol. 191(5), pp. 739-756.

6. Chatelet E., Lornage D. and Jacquet-Richardet G., 2000, Dynamic behaviour of thin-walled composite shafts: A three dimensional approach, 5th annual Engineering System Design and Analysis Conference, ASME, Montreux Switzerland, pp. 1-5

7. Bauchau O., 1983, Optimal Design of High Speed Rotating Graphite/Epoxy Shafts, Journal of Composite Materials, vol. 17(3), pp. 170-181.

8. Gubran H.B.H and Gupta, K., 2005, The effect of stacking sequence and coupling mechanisms on the natural frequencies of composite shafts, Journal of Sound and Vibration, Vol. 282, pp. 231-248.
9. Pereira J. C. and M. E. Silveira, 2002, Evaluation and Optimization of the instability regions on rotors in wounding shaft, II Congresso Nacional de Engenharia Mecânica, João Pessoa.

10. Wettergren H.L. and Olsson K.O., 1996, Dynamic Instability of A Rotating Asymmetric Shaft with Internal Viscous Damping Supported in anisotropic bearings, Journal of Sound and Vibration, Academic Press Limited, Vol. N° 195, pp. 75-84.

11. Newkirk B.L., 1924, Shaft Whipping, General Electric Review, vol. 27(3), pp. 169-178.

12. Kimball A.L., 1925, Internal friction as a cause of shaft whirling, Philosophical Magazine, Series 6, vol. 49(1), pp. 724.

13. Bucciarelli L.L., 1982, On the Instability of Rotating Shafts due to Internal Damping, Journal of Applied Mechanics, Vol. 49, pp. 425.

14. Tsaï, S.W., 1988, Composites Design, 4th edition, Dayton, Ohio, USA.

15. Singh S.P. and Gupta, K., 1994, Damped free vibration of layered composite cylindrical shell, Journal of Sound and Vibration, Vol. 172(2), pp. 191-209.

16. Lalanne M. and Ferraris G., 1998, Rotordynamics Prediction in Engineering, 2nd edition, J. Wiley and Sons

17. R. Sino, T. N. Baranger, E. Chatelet and G. Jacquet, Stability analysis of internally damped rotating composite shaft considering transversal shear, Composites Sciences and Technology, 68, 337-345, 2008.

18. Zinoviev A. and Ermakov N., 1994, Energy Dissipation in Composite Materials, Technic publishing company, Inc., Pennsylvania U.S.A.

19. ROTORINSA®, 2008, "Finite element software for the prediction of the dynamic behavior of rotors in bending", Laboratoire de Mécanique des Contacts et Structures, LaMCoS, UMR5259, INSA Lyon, France.