# Abstract

The physics lying behind rotordynamics is complex to model, so that in many cases numerical processing is the only feasible approach. Being rotordynamics a field of great interest in the aerospace industry, the efforts devoted to its understanding are increasing day by day. Following this tendency, the aim of the present study is to develop a simplified elastodynamic model for the case of rotating structures such that can be addressed through numerical tools, built using the finite element method. For the purpose of analysing the vibration phenomena, modal decomposition and numerical integration have been taken advantage from. In this context, it has been found that the singular value decomposition could be applied in structural analysis to extract dominant displacement fluctuations, allowing the unfolding of global properties of the dynamic response. In the present report, the singular value decomposition has been applied to cantilever beams undergoing a single rotation, giving rise to reasonably satisfactory results.

Keywords: Blades, FEM, rotordynamics, SVD, vibration

# List of Symbols

The next list describes several symbols that will be later used within the body of the document.

#### Abbreviations

• ${\displaystyle {\textrm {CAD}}}$ Computer Aided Design
• ${\displaystyle {\textrm {DOF}}}$ Degree Of Freedom
• ${\displaystyle {\textrm {FEM}}}$ Finite Element Method
• ${\displaystyle {\textrm {SVD}}}$ Singular Value Decomposition
• ${\displaystyle {\textrm {RSV}}}$ Right Side Vector
• ${\displaystyle {\textrm {LSV}}}$ Left Side Vector

#### Mathematical notation

• ${\displaystyle i}$ Imaginary unit
• ${\displaystyle \mathbf {I} }$ Identity matrix
• ${\displaystyle {\mathfrak {L}}}$ Lagrangian
• ${\displaystyle n}$ Outward normal unit vector
• ${\displaystyle \mathbb {R} }$ Set of real numbers
• ${\displaystyle \mho }$ Domain
• ${\displaystyle \Gamma }$ Boundary
• ${\displaystyle \varnothing }$ Null space
• ${\displaystyle \nabla ^{s}}$ Symmetric gradient
• ${\displaystyle {\widetilde {\square }}}$ Spin of a vector
• ${\displaystyle {\square }^{T}}$ Transpose
• ${\displaystyle {\dot {\square }}}$ First temporal derivative
• ${\displaystyle {\ddot {\square }}}$ Second temporal derivative
• ${\displaystyle {\square }^{e}}$ Elemental variable
• ${\displaystyle {\square }_{s}}$Variable expressed in the ${\displaystyle s}$ frame

#### General variables

• ${\displaystyle A}$ Area
• ${\displaystyle \rho }$ Density
• ${\displaystyle \mathrm {F} }$ Force
• ${\displaystyle \mathbf {J} }$ Inertia tensor
• ${\displaystyle m}$ Mass
• ${\displaystyle \mathrm {M} }$ Moment
• ${\displaystyle \mathrm {P} }$ Power
• ${\displaystyle R}$ Radius
• ${\displaystyle t}$ Time
• ${\displaystyle \mathrm {T} }$ Torque

#### Coordinate systems

• ${\displaystyle {\boldsymbol {x}}=\left\{x\,y\,z\right\}^{T}}$ Inertial frame
• ${\displaystyle {\boldsymbol {\varsigma }}=\left\{\varsigma \,\varrho \,\varpi \right\}^{T}}$ Corotational frame
• ${\displaystyle {\boldsymbol {\xi }}=\left\{\xi \,\eta \,\zeta \right\}^{T}}$ Parent domain

#### Kinematics

• ${\displaystyle a}$ Linear acceleration
• ${\displaystyle \alpha }$ Angular acceleration
• ${\displaystyle \theta }$ Angular position
• ${\displaystyle \Omega }$ Angular velocity
• ${\displaystyle x}$ Position
• ${\displaystyle \mathbf {Q} }$ Rotation matrix
• ${\displaystyle \mathbf {T} }$ Rotation tensor
• ${\displaystyle \mathbf {U} }$ Total displacements
• ${\displaystyle \mathrm {v} }$ Velocity

#### Elasticity

• ${\displaystyle k}$ 1D stiffness
• ${\displaystyle \mathbf {u} }$ Elastic displacements in ${\displaystyle {\boldsymbol {\varsigma }}}$
• ${\displaystyle \mathbf {R} }$ Elastic reactions
• ${\displaystyle \mathbf {C} }$ Elasticity tensor
• \mathfrak{u} Generic elastic displacement
• ${\displaystyle \nu }$ Poisson's coefficient
• ${\displaystyle {\overline {\mathbf {t} }}}$ Prescribed tractions
• ${\displaystyle {\overline {\alpha }}\,,\,{\overline {\beta }}}$ Rayleigh damping parameters
• ${\displaystyle \mathbf {W} }$ Reaction torque
• ${\displaystyle I}$ Second moment of area
• ${\displaystyle {\boldsymbol {\varepsilon }}}$ Strain tensor
• ${\displaystyle {\boldsymbol {\sigma }}}$ Stress tensor
• ${\displaystyle E}$ Young's modulus
• ${\displaystyle {\overline {\mathbf {D} }}}$ Viscoelasticity tensor

#### Finite element method

• ${\displaystyle \mathbf {A} }$ Assembly operator
• ${\displaystyle \mathbf {H} }$ Circulatory matrix
• ${\displaystyle \mathbf {CN} }$ Connectivity matrix
• ${\displaystyle \mathbf {D} }$ Damping matrix
• ${\displaystyle \mathbf {J} ^{e}}$ Elemental Jacobian
• ${\displaystyle \mathbf {f} }$ Force vector
• ${\displaystyle \xi _{g}}$ Gauss point
• ${\displaystyle w_{g}}$ Gauss weight
• ${\displaystyle \mathbf {d} }$ Global displacements vector
• ${\displaystyle \mathbf {c} }$ Global variations vector
• ${\displaystyle \mathbf {G} }$ Gyroscopic matrix
• ${\displaystyle \mathbf {M} }$ Mass matrix
• ${\displaystyle \mathbf {N} }$ Matrix of shape functions
• ${\displaystyle \mathbf {B} }$ Matrix of the gradient of ${\displaystyle \mathbf {N} }$
• ${\displaystyle n^{e}}$ Nodes of the element e
• ${\displaystyle n_{el}}$ Number of elements
• ${\displaystyle n_{pt}}$ Number of nodes
• ${\displaystyle n_{d}}$ Problem dimension
• ${\displaystyle {\mathcal {}}{R}}$ Rigid body modes
• ${\displaystyle \mathbf {l} }$ Set of unconstrained DOFs
• ${\displaystyle \mathbf {r} }$ Set of constrained DOFs
• ${\displaystyle \mathbf {K} }$ Stiffness matrix
• ${\displaystyle {\mathcal {}}{V},v_{i}}$ Test functions
• ${\displaystyle {\mathcal {}}{S},u_{i}}$ Trial functions

#### Modal analysis

• ${\displaystyle {\overline {\xi }}_{i}}$ Damping ratio
• ${\displaystyle \lambda _{j}}$ Eigenvalue
• ${\displaystyle \mathbf {q} }$ Generalised coordinates
• ${\displaystyle \omega _{j}}$ Natural frequency
• ${\displaystyle {\boldsymbol {\Phi }}_{j}}$ Natural mode

#### Numeric integration and SVD

• ${\displaystyle \beta }$ Explicit character parameter
• ${\displaystyle \mathbf {U} }$ Matrix of left side vectors
• ${\displaystyle \mathbf {V} }$ Matrix of right side vectors
• ${\displaystyle {\boldsymbol {\Sigma }}}$ Matrix of singular values
• ${\displaystyle \gamma }$ Numerical viscosity
• ${\displaystyle \Delta {\textrm {t}}}$ Timestep

# 1. Introduction

The aim of this section is to briefly introduce the reader to the field of rotordynamics and outline the main models used nowadays to tackle these kind of problems. An overview of the report structuring will also be presented, as well as a review on the actual possibilities to develop the present work. This introduction follows the key concepts presented in [1], [2] and [3] .

## 1.1 Historical frame and background

Since the invention of the potter's wheel, almost 6000 years from now, humanity has relied on rotating machinery for a wide range of activities and processes. This idea of a symmetrical object turning on an axis is the key behind complex inventions such as mills, water wheels and even engines and jet turbines.

Behind those great inventions lies a great understanding of the physics describing this rotational motion. Greeks were the first who tackled machinery in a systematic and even mathematical point of view, bringing new inventions such as Archimedes' screw and Hero's aelopile, the first reaction turbine — which unfortunately could not produce useful work.

 Figure 1.1: Replica of Hero's aelopile [4]

Water wheels became an important source of energy in ancient Chinese civilisation, and later on in the Roman empire as a substitute for manual and animal work. It was used for a wide range of applications, including iron casting and grain grinding. Wind mills took longer to appear, the first practical wind-powered machines dating from the ${\textstyle 9^{th}}$ century in Persia, and reaching Europe three centuries later. These kind of rotating machinery were, however, quite basic, and needed from better scientific background in order to be meaningfully developed.

The world had to wait until the Renaissance Period to see substantial contributions to this field. Leonardo da Vinci is credited for some fundamental contributions to solid mechanics, which were improved by Galileo and later on by Euler and Bernoulli, who formulated the correct equations for simple bending. In the zenith of scientific revolution, Newton's development of calculus unfold a new range of applications of mathematical methods to engineering. It was thanks to the latter that Daniel Bernoulli could formulate the differential equation of motion for a vibrating beam, to be later accepted by Euler when investigating beams under different loading conditions.

The Euler-Bernoulli beam theory was improved by Rayleigh (1877) and later on by Timoshenko (1921), adding the effect of shear. This beam theory, valid even today in the ${\textstyle 21^{st}}$ century, is the backbone of all rotor analysis. The two-dimensional study of structures was carried by Germain and Kirchhoff (1850), and lucid discussions about vibrations of elastic systems were carried by Lord Rayleigh.

Thanks to the contribution of these and many other scientist and mathematicians,the fundamentals of the Theory of Elasticity could be established during the scientific revolution. These deformable bodies were referred as structures by engineers. As the industrial revolution began and the first steam and internal combustion engines appeared, engineers needed from approximate methods, based on energy principles, to carry out accurate studies for solid structures with intricate geometry.

By then, as the required speed of the machinery was rather slow, reciprocating machinery could be used instead of rotating devices. In fact, little was known about the mechanical behaviour of beams under rotating forces. A first study on the topic was carried by Rankine (1869), proposing that a critical speed exists above which rotation becomes impossible. This was later denied by De Laval (1883), a pioneer in the design of steam turbines.

Once the physics of rotordynamics were well understood and steam could be used as motive force, there was a great improvement in the capacity of power generation. The first reaction turbine came up in 1884 at the hands of Charles Parsons. Later on, great development was achieved in the field of rotordynamics due to the invention of the Dynamo in 1878, which led to an outstanding expansion of the steam turbine.

Little time had to pass for rotating machinery to be used as source of propulsion. The first turbojet-powered aircraft was the Heinkel HE-178, dating from 1939. Since then, high-speed complex rotating machinery has taken an important role in many industrial and aeronautical fields. However, they require from precise estimation to work properly, as they are submitted to high thermal and stress loads. Elasticity theories and energy methods are relatively easy to apply to bars, beams, discs, plates and even membranes, but when it comes of complex three-dimensional structures, no analytical solution exists.

In this context, numerical analysis arose. Until the half of the past century, analysis of complex structures was already been carried in a numerical approach using the matrix method. However, even tough it gave good results for simple structures under static conditions, the emergence of aeronautics uncovered some of its main drawbacks, such as the impossibility to describe the flutter phenomena.

 Figure 1.2: Picture of the Heinkel HE-178 [5]

With the release of ${\textstyle NASTRAN}$ by the end of the sixties and further publications pointing its mathematical basis, the finite element method gained popularity among physicists and engineers for its wide range of application. Nowadays, its main drawback is the computational cost, as the computation time increases with the number of elements to simulate.

This problem has already been addressed by applying modal decomposition analysis, which provides great results with less computation cost for dynamics systems under harmonic excitations. Unfortunately, physics describing the vibration phenomena in rotating structures are highly nonlinear, and thus alternative paths are to be sought. The singular value decomposition seems to be a good candidate that, in combination with modal analysis and numerical integration, could unfold predominant patters of the structure. If applied to rotordynamics, some phenomena of aeronautical interest could be studied, including the vibrations generated in helicopters’ and wind turbines’ blades. Other effects caused by the rotation conditions could also be addressed, such as the centrifugal stiffening and resonance phenomena. The latter is, in fact, one of the fields of investigation that lead to the foundation of the finite element method.

## 1.2 State of the art

This section reviews some of the basic topics of rotordynamics and introduces the reader to some essential but important concepts and models related with the dynamics of rotating structures.

### 1.2.1 Introduction to rotordynamics

Rotordynamics is the branch of applied mechanics dealing with devices in which at least one part, the rotor, rotates with significant angular momentum. Although the definition of rotor states that a set of hinges or bearings constrain the position of the rotation axis, there are some cases in which the rotor is not restricted at all. These are known as free rotors, with a spin speed governed by conservation of angular momentum. Examples of free rotors are spinning projectiles, space vehicles and even celestial bodies such as spinning neutron stars.

The aim of this report is, however, to study the so-called fixed rotors, in which the spin axis is somehow fixed and the spin speed imposed by a driving device. First rotordynamic studies date from the late ${\textstyle 19^{th}}$ century, when many methods to deal with rotatory machinery were developed (see A.2). One of the simplest models that described the behaviour of rotatory structures was the Jeffcott rotor, which can qualitatively explain many important features of real-life rotors. However, this simple model fails in computing precise values and in modelling complex machinery, as cannot predict some phenomena such as the dependence of the natural frequencies on the rotational speed.

As time passed by, higher power density, lower weight and faster machinery were demanded. As a consequence, correct quantitative prediction became particularly important. Torque is usually the critical factor when sizing the machine, so in applications related with power generation higher rotation speed is a must. With the increase in speed, power devices became lighter by using stronger but lighter materials. However, stronger doesn't always mean stiffer, an these lighter structures were more prone to vibrate.

Another field that demanded high performing features was propulsion, which needed from high thermodynamic efficiency and power density. That meant that the same amount of heat had to be generated in a smaller volume, and hence with lower thermal capacity. Rotating machines working under these conditions not only have to handle high temperatures, but also stresses due to the thermal gradient.

Although there has been an outstanding progress in the field of rotordynamics during the recent years, it is still a field of very active research. Attitude control of space vehicles, design of wind turbines, industrial rotating machinery, turbojets, electric motors, power generators, helicopter's rotors and space launchers (rockets) are only a few of all the engineering fields that could benefit from those investigations in the near future.

### 1.2.2 Linear rotordynamics

#### Rigid body equations

Imagine a rigid body rotating and moving around, with mass ${\textstyle m}$ and principal moments of inertia ${\textstyle J_{\varsigma },J_{\varrho },{and}J_{\varpi }}$ referred to a frame ${\textstyle \varsigma \varrho \varpi }$ fixed to the body and coincident with its principal axes of inertia. The equations describing its motion are quite complex, with non-linearities in the rotational degrees of freedom. Being ${\textstyle xyz}$ the inertial frame of reference, the equations of motion under the effect of a moment ${\textstyle \mathbf {M} }$ and force ${\textstyle \mathbf {F} }$ are:

 ${\displaystyle \mathbf {F} =\mathrm {m} {\frac {d\;\mathbf {v} }{d\mathrm {t} }}=m\;\mathbf {a} \;\;\;,\;\;\;\mathbf {J} {\dot {\boldsymbol {\Omega }}}+{\boldsymbol {\Omega }}\times (\mathbf {J} {\boldsymbol {\Omega }})=\mathbf {M} }$
(1.1)

Which correspond to Newton's second law and Euler's equation for rigid body dynamics, being ${\textstyle {\boldsymbol {\Omega }}}$ the angular velocity of the body and ${\textstyle \mathbf {a} }$ its linear acceleration. Splitting into components, yields to:

 {\displaystyle \left\{{\begin{aligned}m{\ddot {x}}&=F_{x}\\m{\ddot {y}}&=F_{y}\\m{\ddot {z}}&=F_{z}\end{aligned}}\right.\qquad \left\{{\begin{aligned}M_{\varsigma }&={\dot {\Omega }}_{\varsigma }\,J_{\varsigma }+\Omega _{\varrho }\,\Omega _{\varpi }\,\left(J_{\varpi }-J_{\varrho }\right)\\M_{\varrho }&={\dot {\Omega }}_{\varrho }\,J_{\varrho }+\Omega _{\varsigma }\,\Omega _{\varpi }\,\left(J_{\varsigma }-J_{\varpi }\right)\\M_{\varpi }&={\dot {\Omega }}_{\varpi }\,J_{\varpi }+\Omega _{\varsigma }\,\Omega _{\varrho }\,\left(J_{\varrho }-J_{\varsigma }\right)\end{aligned}}\right.}
(1.2)

As can be spotted, Euler equations are nonlinear in the angular velocity. Moreover, the previous equations only describe the behaviour of a rigid body and elasticity plays no role whatsoever. Thus, one can state that the equations describing the behaviour of an structure under rotational conditions will be highly nonlinear.

However, a linearised model describing the basic features of rotating systems can be obtained embracing some simplifications :

• The rotor has a rotation axis that coincides with one of its baricentrical principal axes of inertia (the rotor is perfectly balanced).
• The deviations from the previous condition are assumed to be small.
• Displacements and velocities (linear and angular) are assumed to be small, with the exception of the rotation angle and angular velocity about the spin axis, which in many cases will be imposed by the driving system.

#### Equation of motion

If the previous simplifications are taken into account, a rotor with constant spin velocity ${\textstyle \Omega }$ can be modelled discreetly by using the adequate elastic model, and the following dynamic equilibrium equation is obtained:

 ${\displaystyle \mathbf {M} \,{\ddot {\mathbf {q} }}(\mathbf {t} )+(\mathbf {D} +\mathbf {G} )\,{\dot {\mathbf {q} }}(t)+(\mathbf {K} +\mathbf {H} )\,\mathbf {q} (t)=\mathbf {f} (t)}$
(1.3)

In the previous equation, some of the present terms may look familiar to the reader: ${\textstyle \mathbf {M} }$, ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {D} }$ stand for the symmetric mass, stiffness and damping matrices, while ${\textstyle \mathbf {q(t)} }$ is a vector containing the generalised coordinates referred to the inertial frame ${\textstyle xyz}$ and ${\textstyle \mathbf {f(t)} }$ a time-dependant vector that accounts for the effect of external forces. However, in the context of rotordynamics, new matrices appear: ${\textstyle \mathbf {G} }$ and ${\textstyle \mathbf {H} }$ are the skew-symmetric gyroscopic and circulatory matrices, respectively.

The gyroscopic matrix contains inertial and conservative terms, linked with the gyroscopic moments that act on the rotating elements of the structure. If the equations are referred to the rotating and thus noninertial frame of reference, the Coriolis effect is included in the gyroscopic matrix. The circulatory matrix, on the other hand, contains the nonconservative effects related with internal damping of the rotating parts and, in some cases, accounts for the effect of the surrounding fluids. In many cases, it is found that in the field of rotordynamics, both ${\textstyle \mathbf {G} }$ and ${\textstyle \mathbf {H} }$ are proportional to the spin velocity ${\textstyle \Omega }$, and when the latter tends to zero, the gyroscopic and circulatory skew-symmetric elements vanish and the system reduces to that of a static structure.

It is important to notice that in this context, ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {D} }$ may not have the same expression they had in the case of still strictures. In fact, they can be split into two matrices: one corresponding to static structure and another which depends on the rotating speed, often on ${\textstyle \Omega ^{2}}$.

#### Complex coordinates

In some cases, the simplifications made in the previous page cannot be assumed, and the study of the rotor becomes very complicated. In this situations, however, an model resembling equation 1.3 can be obtained, although with reference to a noninertial frame. Moreover, most flexible rotors can be assumed to behave as beams. With these hypotheses, the lateral behaviour of the rotor can be uncoupled from its axial and torsional, simplifying thereby the resulting equations. Mathematically, this is done by using complex coordinates.

If the problem can be reduced into two-dimensional, that is, spinning around the ${\textstyle z}$ axis and displacements only in the ${\textstyle xy}$ plane, the latter can be expressed in terms of a complex number ${\textstyle r(t)=x(t)+iy(t)}$, where ${\textstyle i}$ stands for the imaginary unit ${\textstyle {\sqrt {-1}}}$. This formulation is equivalent to representing displacements in vectorial form, but allows a more convenient analytical form to solve the rotor equation of motion (eq. 1.3).

The key of this method is that symmetric matrices are real when the equations are written in terms of complex coordinates, whereas skew-symmetric matrices become symmetric but with imaginary terms. Using the complex number formulation, the equation of motion for a axisymmetrical rotor results in:

 ${\displaystyle \mathbf {M} ^{\prime }\,{\ddot {\mathbf {q} }}^{\prime }(t)+\left(\mathbf {D} ^{\prime }+i\mathbf {G} ^{\prime }\right)\,{\dot {\mathbf {q} }}^{\prime }(t)+\left(\mathbf {K} ^{\prime }+i\mathbf {H} ^{\prime }\right)\,\mathbf {q} ^{\prime }(t)=\mathbf {f} ^{\prime }(t)}$
(1.4)

Separating real from imaginary parts:

 ${\displaystyle \left[{\begin{array}{cc}{\mathbf {M} ^{\prime }}&{\mathbf {0} }\\{\mathbf {0} }&{\mathbf {M} ^{\prime }}\end{array}}\right]{\ddot {\mathbf {q} }}(t)+\left(\left[{\begin{array}{cc}{\mathbf {D} ^{\prime }}&{\mathbf {0} }\\{\mathbf {0} }&{\mathbf {D} ^{\prime }}\end{array}}\right]+\left[{\begin{array}{cc}{\mathbf {0} }&{\mathbf {-G'} }\\{\mathbf {G'} }&{\mathbf {0} }\end{array}}\right]\right){\dot {\mathbf {q} }}(t)+\left(\left[{\begin{array}{cc}{\mathbf {K} ^{\prime }}&{\mathbf {0} }\\{\mathbf {0} }&{\mathbf {K} ^{\prime }}\end{array}}\right]+\left[{\begin{array}{cc}{\mathbf {0} }&{\mathbf {-H'} }\\{\mathbf {H'} }&{\mathbf {0} }\end{array}}\right]\right){\mathbf {q} }(t)={\mathbf {f} }(t)}$
(1.5)

The solution of the previous system can be written in terms of the complex frequency ${\textstyle \;\;\mathbf {s=\sigma {+}i\omega } }$: the imaginary part ${\textstyle \omega }$ stands for the frequency of the free motion (whirl frequency) while the real part ${\textstyle \sigma }$ is the decay rate changed in sign.

The main advantage of this method is that allows the system to be analysed using modal decomposition, which in conventional coordinates required from all the matrices to be square and symmetrical. Further details on the modal analysis and computation of eigenvectors and eigenvalues can be found in section 5.1.

#### The Campbell diagram

As has been seen, the spin speed of the structure can appear explicitly in the equation of motion. Thus, a dependence of the natural frequencies on this speed is expected. If that is the case, it is quite common to plot the natural frequencies ${\textstyle \omega }$ as functions of ${\textstyle \Omega }$. In many cases, frequencies characterising exciting forces also depend on ${\textstyle \Omega }$ and so are reported in the same graph, which is commonly known as Campbell diagram. Clearly, it is symmetrical for both ${\textstyle \omega }$ and ${\textstyle \Omega }$ axis, so only the first quadrant is represented. The intersections with the ${\textstyle \omega }$ axis correspond to the natural frequencies at standstill of the system (${\textstyle \Omega }$ = 0). It is important to note that the Campbell diagram can only be graphed if the equation of motion is linearised, because only in this case the concept of natural frequencies applies.

If the exciting forces are time-dependent with an harmonic time story, they can be plotted in the Campbell diagram. For instance, that is the case of the forces resulting from the unbalance of the rotor, which can be modelled with frequency equal to the rotation speed ${\textstyle \Omega }$. If the shape of the force function is less regular, it can be split into harmonic terms using Fourier decomposition. The complete response of the system will be the sum of the responses of the system under those harmonic components.

Quite often, the relation between the frequency of the forcing function and ${\textstyle \Omega }$ is of simple proportionality, and thus can be represented in the Campbell diagram by a straight line. If the proportion is one to one, the line ${\textstyle \omega =\Omega }$ is graphed, and the excitation is said to be synchronous. The spin speeds ${\textstyle \Omega }$ at which a forcing function has a frequency coinciding with a natural frequency of the system at that same ${\textstyle \Omega }$ are known as critical speeds. They can be found in the diagram by looking for the intersections between natural and forcing frequencies.

As one may suspect, the critical speeds bring up resonance, and thus may be avoided as great displacements may be encountered. If that is the case, the rotor cannot operate near this speed without encountering strong vibrations and even catastrophic failure. However, not all critical speeds are equally dangerous, and that will depend on the main modes of the structure. Flexural critical speed are those linked with the flexural natural frequencies, and use to be particularly dangerous. They can be detected in the Campbell diagram by looking for the intersection with the line ${\textstyle \omega =\Omega }$. In addition to these flexural speeds, torsional critical speeds can also be very dangerous.

The range starting from zero to the first critical speed is referred as subcritical range, and above it, supercritical range starts. In the past, it was thought that operation above the first critical speed was impossible, but it was later denied by De Laval, Stodola and many other rotordynamists. It is easy to mix up natural and critical speeds, mostly due that in many cases their value coincides. Even if that happens, the two physical phenomena are very different, particularly concerning the stressing of the rotor.

 Figure 1.3: Example of Campbell diagram [6]

To illustrate a little bit this section, Figure 1.3 shows the Campbell diagram of a compressor rotor blade. Bold curves represent natural frequencies of the structure, while the straight lines show the exciting frequencies. Each intersection between them is a resonance point. Here, rotation introduces an stiffening effect: natural frequencies increase with the spin speed.

It is important to remember that the concept of critical speeds and the Campbell diagram only apply to linear systems, although a more arbitrary definition of these speeds can be used in the case of nonlinear systems, referring the speed at which greater amplitudes are found.

#### Frequency and time domains

Linear rotordynamics can be analysed in the frequency domain, a widely used approach that brings up an insight on how the rotor behaves under different conditions. These kind of solutions allow the engineers to perform parametric design and optimization studies, as natural frequencies and vibration modes are easily obtained using this frame of analysis. When the system is linear and operating in steady-state condition, the most popular tool to analyse its behaviour is modal decomposition, which reduces the system of equations to an eigenvalue problem.

This is achieved, however, when the rotating structure is idealised. If the system is studied in the transient state, no frequency analysis can be achieved and only numerical integration in time is possible. Although integrating numerically the equations of motion is a widely used method nowadays to tackle complex problems, it only allows to solve very particular cases without bringing up an overall insight on the phenomena. Moreover, there is a high computational cost associated to this kind of numerical methods which is not present when performing a frequency analysis.

### 1.2.3 Nonlinear and nonstationary rotordynamics

The previous section presented a linearised model for a rotating structure. This is, however, only an idealisation as real rotors always deviate from the lineal model. Rotors are often designed to behave linearly in their nominal conditions. However, some sources of nonlinearity such as bearings, dampers and seals may disturb the system from its ideal behaviour. Typical approaches are to consider those elements as rigid bodies, or to linearise its behaviour, in order to simplify the resulting equations. This, however, leads to unsatisfactory approximations.

If the structure is highly nonlinear, sub- and super-harmonics play an important role, as well as chaotic oscillations. Chaos has been spotted in nonlinear rotor models, and numerical results indicate that they are present in turbojet engines. If the whirl amplitude grows in time, nonlinearities emerge until the system reaches a limit cycle.

As is to expect from nonlinear systems, sudden jumps from an equilibrium configuration to another may happen. Moreover, hysteresis may appear, making the spin-up of the system different from the deceleration, even if the speeds are the same. When the system to study becomes that complex, the only way to solve the equation is by means of numerical integration.

The previous linear analysis required from constant spin velocity ${\textstyle \Omega }$ to be valid. It is, in fact, a common simplification in rotordynamics analysis. This model do not explains, however, the accelerating and decelerating transient states, or quick changes in the forcing functions. Even if the rotor is assumed to be linear, angular accelerations may introduce new terms in the matrices of equation 1.3 which will also become time-dependant. As a consequence, the system becomes complicated and the equations of motion are no longer solvable analytically. As an integrating numerical scheme will be required to solve the nonstationary system, there is little advantage in previously linearising the equation of motion.

### 1.2.4 The Jeffcott rotor

The analytical study of rotors dates back to the late nineteenth century, when Föppl, Stodola and Belluzzo described the behaviour of a simplified rotor. This model was deeply studied by Jeffcott (1919), and thus received its name. The model consists of a point mass attached to a massless shaft, being one of the simplest models used to describe the flexural behaviour of rotors. Yet, this is an oversimplification and one has to understand its limitations, although it may bring up qualitative insight into important phenomena of rotordynamics.

The model covered in this section considers an axisymmetrical fixed-free rotor, in which damping plays no role whatsoever. The restoring force is provided by an overall stiffness ${\displaystyle k}$, which will be considered the stiffness of the shaft. Let's consider an inertial frame ${\textstyle xy}$ and rotating noninertial frame ${\textstyle \varsigma \varrho }$. In addition, small displacements are assumed, and the point mass is always to be contained in the ${\displaystyle xy}$ plane. As in classic rotordynamics, the spin speed ${\textstyle \Omega }$ is assumed to be constant.

Previously to the study of the Jeffcott rotor, let's consider an even simpler model consisting in a point mass constrained to move only in the axial direction ${\textstyle \varsigma }$. This is the kind of model Rankine used to develop his rotor analysis. The Lagrangian to the system yields to:

 ${\displaystyle {\mathfrak {L}}={\frac {1}{2}}m\left({\dot {\varsigma }}^{2}+\Omega ^{2}\varsigma ^{2}\right)-{\frac {1}{2}}k\varsigma ^{2}}$
(1.6)

And the equation of motion results in

 ${\displaystyle m{\ddot {\varsigma }}+\left(k-m\Omega ^{2}\right)\varsigma =0}$
(1.7)

This equations shows a decrease in the stiffness known as centrifugal softening, as the term multiplying the displacement is ${\textstyle k-m\Omega ^{2}}$, which is clearly smaller than ${\textstyle k}$. Moreover, this system becomes unstable above the critical speed ${\textstyle \omega ={\sqrt {(k/m)}}}$, prohibiting operation in the supercritical range.

 Figure 1.4: Rankine rotor model [3]

If the guide is removed, displacement in the transversal direction ${\textstyle \varrho }$ becomes possible, and a Jeffcott rotor is obtained. The equation of motion in the rotating coordinates results in

 ${\displaystyle m\left[{\begin{array}{cc}{1}&{0}\\{0}&{1}\end{array}}\right]\left\{{\begin{array}{l}{\ddot {\varsigma }}\\{\ddot {\varrho }}\end{array}}\right\}-2m\Omega \left[{\begin{array}{cc}{0}&{1}\\{-1}&{0}\end{array}}\right]\left\{{\begin{array}{c}{\dot {\varsigma }}\\{\dot {\varrho }}\end{array}}\right\}+\left[{\begin{array}{cc}{k-m\Omega ^{2}}&{0}\\{0}&{k-m\Omega ^{2}}\end{array}}\right]\left\{{\begin{array}{l}{\varsigma }\\{\varrho }\end{array}}\right\}=0}$
(1.8)

This model now takes into consideration Coriolis forces, which push the mass out of is axis, allowing self-centring and thus operation in the supercritical range. The equation of motion can be written in the inertial frame using a simple change of coordinates consisting in a rotation of angle ${\textstyle \Omega t}$. The equation of motion of the Jeffcott rotor in the inertial frame reads as

 ${\displaystyle {\begin{array}{l}{m{\ddot {x}}+kx=0}\\{m{\ddot {y}}+ky=0}\end{array}}}$
(1.9)

The equation above predicts no softening effect nor instability in the inertial frame, and a whirl frequency ${\textstyle \omega _{n}=\pm {\sqrt {k/m}}}$, equal to the flexural critical speed ${\textstyle \Omega _{cr}}$. For the case of a simply supported shaft, transverse stiffness ${\textstyle k}$ is equal to ${\textstyle {\frac {48EI}{L^{3}}}}$ [7], where ${\textstyle E}$ stands for the Young’s modulus, ${\textstyle I}$ for the second moment of area of the shaft cross-section and ${\textstyle L}$ for the span of the shaft. The above equation predicts that, if the initial conditions are such that the mass moves into circular orbits, the system is not affected by any kind of vibration but describes whirling motion instead. In the inertial frame, angular frequency may not be constant, but the period will be equal to ${\textstyle 2\pi /\omega _{n}}$.

The natural frequency of the Jeffcott rotor does not depend on ${\textstyle \Omega }$, and thus the Campbell diagram consists in horizontal straight lines. This is the simplest way of tackling the Jeffcott rotor, although more realistic variations of the model are widely used, taking into account forcing functions, damping, shaft bow, unbalance and out-of-plane displacements.

### 1.2.5 Gyroscopic effect

In the previous section, the Jeffcott model considered a rotor consisting in a point mass, and consequently, its moments of inertia were not taken into account. This is, however, a huge simplification. In reality, Campbell diagrams are not horizontal lines, but natural frequencies of bending modes depend on the spin speed. If one aims to account for this behaviour, has to consider the influence of inertia moments. Strictly speaking, the system has six degrees of freedom, and so six generalised coordinates should be defined in order to study its dynamic behaviour.

This behaviour was first noticed by Rayleigh, who identified the effect of the rotatory inertia of a disk mounted on a shaft. Stodola also considered this effect, and identified it as the responsible of giving rise to the split of natural frequencies. The study of critical speeds under the gyroscopic effect was carried by Green, and later studied using energy methods by Carnegie.

One of the simplest models follows the main assumptions of the Jeffcott rotor, but with a rigid body with nonvanishing moments of inertia instead of a point mass. In the undeformed position, one of the principal axes of inertia coincides with ${\textstyle z}$. Principals moments of inertia will be referred as polar moment of inertia ${\textstyle J_{p}}$, about the rotation axis, and transversal moment of inertia ${\textstyle J_{t}}$, about any axes in the rotation plane. In matrix form:

 ${\displaystyle \mathbf {J} =\left[{\begin{array}{ccc}{J_{t}}&{0}&{0}\\{0}&{J_{t}}&{0}\\{0}&{0}&{J_{p}}\end{array}}\right]}$
(1.10)
 Figure 1.5: Massless shaft modes (lumped mass and disk models)

In order to illustrate the problem, imagine a beam and its two first structural modes (Figure 1.5). The first mode introduces a kinetic energy term, as the mass and the disk have been displaced by the deflection of the shaft. However, if the shaft is massless, there is no kinetic energy associated to the second mode, as the point mass has not changed its position. On the contrary, if the disk rotation given by the slope ${\textstyle (\partial v/\partial x)}$ is taken into account, kinetic energy associated to the angular velocity emerges and has the following expression:

 ${\displaystyle {\frac {1}{2}}\;J_{T}\left\{{\frac {d}{dt}}\left({\frac {\partial v}{\partial x}}\right)\right\}^{2}}$
(1.11)

The above value of rotation energy is clearly zero if the beam is considered to be a point mass with no inertia, but when dealing with disks and rigid bodies, the inertia term cannot be neglected and neither can the rotation energy associated. This was a very known problem in static structures, known as rotatory inertia of beams. However, when the shaft rotates, the body acts like a spin top and gyroscope, having a significant effect in rotordynamics.

To understand the physics behind this phenomena, one has to remember Euler's equation for rigid bodies, which are clearly nonlinear with the angular velocities (see equation 1.2). Imagine a freely spinning disk (for instance, Figure 1.6) with angular velocity ${\textstyle \Omega }$, rotating around an axis perpendicular to its non-deformed symmetry plane. Consider now that a whirling motion is introduced now around the ${\textstyle x}$ axis, due to any kind of perturbation or imbalance.

 Figure 1.6: Free spinning disk [1]

Because of whirling, precessional motions linked to the shaft slopes at the disk centre are introduced. The disk then moves similar to a spin top, and the gyroscopic effect is introduced. Moreover, precession can happen in both ${\textstyle xz}$ and ${\textstyle yz}$ planes, leading to two different gyroscopic torques at which inertia torques (${\textstyle J_{T}\cdot {\ddot {\theta }}_{i}}$ and ${\textstyle J_{T}\cdot {\ddot {\phi }}_{i}}$) are added. In this context, the moment relations across the disc ${\textstyle i}$ are written as

 ${\displaystyle {\begin{array}{l}{M_{y\;i}=J_{p}\cdot \Omega \cdot {\dot {\phi }}_{i}+J_{T}\cdot {\ddot {\theta }}_{i}}\\{M_{z\;i}=-J_{p}\cdot \Omega \cdot {\dot {\theta }}_{i}+J_{T}\cdot {\ddot {\phi }}_{i}}\end{array}}}$
(1.12)

This relations have to be taken into account when modelling the dynamics of the system. In the model seen in equation 1.3, these effects are included in the gyroscopic matrix ${\textstyle \mathbf {G} }$. A further simplification can be achieved, uncoupling rotational and translational motions, leading only to four degrees of freedom. If that is the case, one can define the non-dimensional term ${\textstyle \delta ={\frac {J_{p}}{J_{t}}}}$, and study the behaviour of the rotor as a function of ${\textstyle \delta }$.

The Jeffcott model, for instance, corresponds to ${\textstyle \delta =0}$, in which case curves in the Campbell diagram are horizontal lines. If ${\textstyle \delta >1}$, the rotor is disk-like, and as there is no intersection of the curves in the Campbell diagram with the line ${\textstyle \omega =\Omega }$, no critical speed exist. On the contrary, if ${\textstyle \delta <1}$, the rotor is considered to be very long, and a critical speed linked with conical motion exists. Finally, if ${\textstyle \delta }$ is exactly one, the rotor is spherical and although there is no intersection with the line ${\textstyle \omega =\Omega }$, it tends asymptotically to it as ${\textstyle \Omega }$ increases. Regarding the frequency of whirling, it decreases in absolute value as ${\textstyle \Omega }$ increases for the case of backward whirling, and increases with an inclined asymptote ${\textstyle \omega =\delta \Omega }$ in the case of forward whirling.

### 1.2.6 Centrifugal softening and stiffening

In previous sections, the concepts of centrifugal softening and stiffening have been used to refer the decrease and increase, respectively, of the natural frequencies of a rotating structure, compared with the natural frequencies of the same structure but standstill. While centrifugal stiffening is a very well-known effect, centrifugal softening is "an alleged and elusive phenomena, and some doubts can be cast on the mathematical models causing it to appear" [3]. A model of this type is used for instance in [1] (see page 285), trying to account for the softening effect in a solid rotor. The resulting equation is:

 ${\displaystyle \mathbf {M} \,{\ddot {\mathbf {q} }}+(\mathbf {K} -\Omega ^{2}\mathbf {M} )\,{\mathbf {q} }=0}$
(1.13)

There is a term proportional to the mass matrix and to the square of the spin velocity, which is subtracted from ${\textstyle \mathbf {K} }$ and thus natural frequencies wane. However, when flexibility and inertia of the disks and blades are accounted for, centrifugal stiffening plays a crucial role. It can be explained as the result of the stressing caused by tensile forces, which produces a virtual increase of stiffness proportional to the centrifugal force and thus to ${\textstyle \Omega ^{2}}$. The simplest way to model a system accounting for both softening and stiffening effects is to use a 1${\textstyle \mathbf {\frac {1}{2}} }$ dimensional approach, which consists in modelling shafts as beams and using axisymmetrical elements. The linearised homogeneous equation of motion [3] for the free vibration, taking into account axial, torsional, out-of-plane and in-plane bending behaviour reads as

 ${\displaystyle \mathbf {M\,{\ddot {q}}} +\Omega \,\mathbf {G} \,{\dot {\mathbf {q} }}+\left(\mathbf {K} +\mathbf {K} _{\Omega }\,\Omega ^{2}-\mathbf {M} _{n}\,\Omega ^{2}\right)\,\mathbf {q} =\mathbf {0} }$
(1.14)

${\textstyle \mathbf {M} }$, ${\textstyle \mathbf {G} }$ and ${\textstyle \mathbf {K} }$ are the mass, gyroscopic and stiffness matrices. ${\textstyle \mathbf {K_{\Omega }} }$ is the centrifugal stiffening matrix, a geometric matrix whose effects are proportional to ${\textstyle \Omega ^{2}}$. On the other hand, ${\textstyle \mathbf {M_{n}} }$ is the centrifugal softening matrix, responsible for causing a decrease of the natural frequencies. The key is to assess whether natural frequencies actually decrease or not as ${\textstyle \Omega }$ increases. To check that, it must be known whether the effect of ${\textstyle \mathbf {M_{n}} }$ is stronger than ${\textstyle \mathbf {K_{\Omega }} }$ or not.

The effect of ${\textstyle \mathbf {K_{\Omega }} }$-${\textstyle \mathbf {M_{n}} }$ is dual: a stiffening effect due to tensile stresses caused by rotation and a softening due to the effects of the noninertial frame in which the equations are written. It can be demonstrated (see [3]) that that matrix is positive defined, and thus the overall net effect is of stiffening. The only models that, for its simplicity, are not able to capture the stiffening effect are the Rankine and Jeffcott rotors, which neither account for the gyroscopic effect. Moreover, the Rankine model is the only one that predicts a strong softening effect, which is known to lead to incorrect results.

## 1.3 Report outline

As presented in the previous section, the dynamic behaviour of rotormachinery is complex and the reigning equations can easily become convoluted, leading to a set of equations whose solving is far out of the scope of this report. Hence, instead of modelling the structure with complex and accurate equations, simpler models will be considered. The goal is not to get credible results, but to study on the methods used to solve the equations. For the cases where ${\textstyle \Omega }$ is not constant, the equations cannot be solved in terms of modal decomposition and only numerical integration is possible. These numerical schemes take a lot of time and computation cost, and give only particular solutions of the system. The key of this approach is to develop a numeric method based on the SVD that has the advantage of the modal analysis of capturing global properties of the system —such as vibration modes and natural frequencies— when applied to rotors with angular acceleration and other nonlinearities.

Regarding the written organization of the document, it is worth noting that an appendix is included at the end of the report (sections enumerated with capital Latin letters). Particular topics that could be from interest to the reader are presented there, as well as some figures related with the simulations carried out in chapters 6 and 7.

Focusing on the report, the chapters into which is divided are:

Chapter 2. The elastic model that describes the vibration behaviour of rotors is formulated. The strong form of the elastodynamic equation is obtained and analytically solved for the 1D stationary case.

Chapter 3. The weak form of the elastodynamic equations is obtained, posed in terms of the finite element method and formulated in the elemental frame. General concepts regarding the FEM methods are reviewed in this section.

Chapter 4. The general steps of every FEM simulation — preprocessing, solver and postprocessing — are reviewed. The 1D stationary case is solved using the FEM.

Chapter 5. Different methods to tackle the dynamic case are compared: modal decomposition analysis, numerical integration and singular value decomposition.

Chapter 6. The methods reviewed in chapter 5 are used to perform the structural analysis of rotating beams. In particular, the chapter focuses on the application of the SVD to recover predominant information of rotating structures.

Chapter 7. The dynamic model is further generalised, including aerodynamic forces and rigid-elastic coupling. The numeric solver is used to simulate the starting of an helicopter rotor.

Chapter 8. The developed work is reviewed and final conclusions are stated.

The presented outline may seem confusing for those with little experience in numerical methods. If that is the case, the reader should keep in mind the following scheme that encompasses the different parts of the numerical tool that is to be developed, as well as its chronological order (Figure 1.7).

 Figure 1.7: Overall numeric scheme 2

(1)

(2) The used software is property of the following trademarks: ${\displaystyle \mathrm {3D} }$ ${\displaystyle \mathrm {SOLIDWORKS} ^{\mbox{™}}}$. 2019. Version 2018 SP4. Dassault Systèmes ${\displaystyle ^{\mbox{®}}}$, Vélizy-Villacoublay, France.

${\displaystyle \mathrm {GiD} ^{\mbox{®}}}$. 2019. Official version 14.0.2. CIMNE ${\displaystyle ^{\mbox{©}}}$, Barcelona, Spain.

${\displaystyle MATLAB^{\mbox{®}}}$. 2019. Version R2018a. The MathWorks, Inc., Natick, MA, USA.

# 2. Rotor modelling

The first step in the study of vibrations in rotating structures is to define an adequate model that describes its behaviour. In section 1.3, a justification on why a simple model is preferable to more accurate but complex equations is presented. The aim of this section is to keep going with this justification, presenting the chosen model and its limitations. The equations of motion will also be presented latter in this chapter.

## 2.1 Problem formulation and hypotheses

As it has already been stated in section 1.2, the study of rotordynamics can easily become complicated, as the understanding of the whole phenomena of rotating structures requires from a deep knowledge in several fields. However, an accurate study of rotatory structures is far from the scope of this project. Instead, more emphasis will be given to the methods used to solve the resulting equations, rather than on the modelling and simulation of the rotor. That is why the model will be defined under the following conditions and hypotheses:

• The axis of rotation of the rotor is known and fixed. Angular velocity and acceleration are known and given by the driving device.
• The problem is restricted to small displacements but great rotations induced by the driving device. Hence, it can be modelled using the linear elastic equations (A.1).
• The rotor will be modelled in terms of the finite element method using solid elements, being nodal displacements the degrees of freedom.
• The inertia of the rotor will not be taken into consideration, and thus no gyroscopic effect is expected.

The first point gives an important hint about the nature of the problem. In the inertial frame, rotating velocities and accelerations will be known. In other words, the inertial kinematic variables are not the unknown, which in fact is the torque needed to deliver the required power to maintain the given angular conditions. One has to take into account that, although the driving engine may introduce large rotations, strains will still be small in a rotating frame of reference static with respect the rotor.

The strongest limitation of the model is given by the forth point: no gyroscopic effect will be taken into consideration. This limits the correctness of the model when applied to structures with non-negligible inertia. The gyroscopic effect would induce rotations in different axes, breaking the first hypothesis. With this type of formulation, great part of the phenomena description is lost, as complex models such as the one of equation A.2 are difficult to model using classical formulations.

The above hypotheses will make the model inherently inaccurate. The resulting simulations will not describe the vibration behaviour of real rotor blades, although it is a great first step into the actual rotordynamic problem. The key of this approach is to make equations easy to pose, and thus expedite the programming of the model in terms of the FEM. Hopefully enough, this will allow a deeper study on the methods used to solve the dynamic system.

## 2.2 Frames of reference and rotation matrix

Previous to the posing of the equations describing the kinematics of the system, it is convenient to define the different frames of references that will be used: an inertial frame ${\textstyle xyz}$, fixed to a static point ${\textstyle O}$; and a noninertial rotating frame ${\textstyle \varsigma \varrho \varpi }$ fixed to the rotor — thus sharing its angular velocity ${\textstyle \Omega }$ and acceleration ${\textstyle \alpha }$. The planes defined by ${\textstyle xy}$ and ${\textstyle \varsigma \varrho }$ are meant to be coincident, and thus ${\textstyle z}$ and ${\textstyle \varpi }$ will be the same axis (since both basis are positive definite). From now on, they both will be referred as ${\textstyle z}$, the rotation axis.

 Figure 2.1: Vector rotation in the plane xy

To change from one coordinate system to another, one can define a rotation matrix, taking advantage of the fact that the angular displacement is known for every value of time. Following Euler's definition of rotation, it is stated that rotation in space can always be described by a rotation along a certain axis over a certain angle. In this case, the axis of rotation is ${\textstyle z}$. One can, for instance, represent the rotation in a two-dimensional space, being ${\textstyle {\mathbf {r} }_{i}}$ and ${\textstyle {\mathbf {r} }_{f}}$ the vectors describing initial and final positions, respectively, of a given point. The rotation is represented in Figure 2.1:

 ${\displaystyle {\mathbf {r} _{f}}={\begin{pmatrix}x_{f}\\y_{f}\end{pmatrix}}={\begin{pmatrix}R\cos(\theta +\varphi )\\R\sin(\theta +\varphi )\end{pmatrix}}={\begin{pmatrix}{\overset {x_{i}}{\overbrace {R\cos(\varphi )} }}\cos(\theta )-{\overset {y_{i}}{\overbrace {R\sin(\varphi )} }}\sin(\theta )\\{\underset {y_{i}}{\underbrace {R\sin(\varphi )} }}\cos(\theta )-{\underset {x_{i}}{\underbrace {R\cos(\varphi )} }}\sin(\theta )\end{pmatrix}}={\begin{bmatrix}\cos(\theta )&-\sin(\theta )\\\sin(\theta )&\cos(\theta )\end{bmatrix}}{\begin{bmatrix}x_{i}\\y_{i}\end{bmatrix}}={\mathbf {{Q}\;_{(\theta )}} }\;{\mathbf {r} _{i}}}$
(2.1)

Where ${\textstyle {\mathbf {{Q}\;_{(\theta )}} }}$ stands for the rotation matrix, and it is satisfied that ${\textstyle \left\|{\mathbf {r} }_{i}\right\|=\left\|{\mathbf {r} }_{f}\right\|=R}$.

It will be interesting, when posing the kinetic equations, to consider the cross product or spin of a vector ${\textstyle t}$, represented with ${\textstyle {\widetilde {t}}}$ and defined as the linear operator:

 ${\displaystyle {t}={\begin{pmatrix}t_{1}\\t_{2}\\t_{3}\end{pmatrix}}\;\Rightarrow \;{\widetilde {t}}={\begin{bmatrix}0&-t_{3}&t_{2}\\t_{3}&0&-t_{1}\\-t_{2}&t_{1}&0\end{bmatrix}}}$
(2.2)

Using this operator, the rotation can be expressed in terms of the Rodrigues formula [8]:

 ${\displaystyle {\mathbf {{Q}\;_{(\theta ,{n})}} }=\left[\cos {\theta }\,{\hbox{I}}+(1-\cos {\theta }){\hbox{n}}{\hbox{n}}^{T}+\sin {\theta }\;{\widetilde {\hbox{n}}}\right]}$
(2.3)

Which for the case ${\textstyle {\mathbf {n} }={\widehat {\mathbf {k} }}}$ gives the same result as in equation 2.1. What is interesting about the spin operator is the easiness to represent matrices derivatives. Defining the temporal derivative of a function as ${\textstyle {\frac {\mathrm {d} \mathbf {f} }{\mathrm {d} t}}={\dot {\mathbf {f} }}}$, one can compute the temporal derivative of the three-dimensional rotation matrix using the spin operator:

 ${\displaystyle {\frac {\mathrm {d} {\mathbf {{Q}\;_{(\theta )}} }}{\mathrm {d} t}}={\mathbf {{\dot {Q}}\;_{(\theta )}} }={\begin{bmatrix}-\sin(\theta ){\dot {\theta }}&-\cos(\theta ){\dot {\theta }}&0\\\cos(\theta ){\dot {\theta }}&-\sin(\theta ){\dot {\theta }}&0\\0&0&0\end{bmatrix}}={\mathbf {{Q}\;_{(\theta )}} }{\begin{bmatrix}0&-{\dot {\theta }}&0\\{\dot {\theta }}&0&0\\0&0&0\end{bmatrix}}={\mathbf {{Q}\;_{(\theta )}} }\;{\widetilde {\dot {\theta }}}\;=\;{\mathbf {{Q}\;_{(\theta )}} }\;\;{\widetilde {\Omega }}}$
(2.4)

The above expression can be generalised to compute matrix derivatives of order ${\textstyle n}$:

 ${\displaystyle {\frac {\mathrm {d^{n}} {\mathbf {{Q}\;_{(\theta )}} }}{\mathrm {d} t^{n}}}={\mathbf {{Q}\;_{(\theta )}} }\;\;{\widetilde {\Omega }}\;^{n}}$
(2.5)

## 2.3 Kinematic equations

The kinematic description of the rotor is going to be quite simple. Imagine a particle of the rotor, represented by a point mass with no inertia associated. This point mass is subjected to a rotation given by the driving device. If the body were perfectly rigid, displacements in the rotating frame ${\textstyle \varsigma \varrho }$ would be zero, and positions in the ${\textstyle xy}$ plane would be given by the driving device. This is, however, not the case. Displacements in all ${\textstyle \varsigma }$, ${\textstyle \varrho }$ and even ${\textstyle z}$ directions may occur.

To illustrate this effect, imagine a two-dimensional case in which elastic displacements are only possible in ${\textstyle \varsigma }$. Imagine the mass point starting from an initial position ${\textstyle {\mathbf {r} }_{i}}$, which after a certain time holds a generic position ${\textstyle {\mathbf {r} }}$. The total displacement ${\textstyle {\mathbf {U} }}$ can be split in terms of rigid body ${\textstyle {\mathbf {U} }\;_{rb}}$ and elastic ${\textstyle {\mathbf {U} }\;_{el}}$ displacements:

 Figure 2.2: Displacement as composition of rotation and deformation

The equations describing this motion are easy to deduce. For simplicity, given a vector ${\textstyle v}$, let's denote as ${\textstyle v_{(x)}}$ the vector described in inertial coordinates, and ${\textstyle v_{\mathbb {(\varsigma )} }}$ the same vector but in the corotational coordinates. For instance, ${\textstyle {\mathbf {r} }_{\;(x)}=\left(x\,y\,z\right)^{T}}$ and ${\textstyle {\mathbf {r} }_{\;(\varsigma )}=\left(\varsigma \,\varrho \,\varpi \right)^{T}}$. It can be easily deduced from equation 2.2:

 ${\displaystyle {\mathbf {U} }\;_{(x)}={\mathbf {r} }_{\;(x)}-{\mathbf {r} }_{\;i\;(x)}={\mathbf {U} }\;_{rb\;(x)}\;+\;{\mathbf {U} }\;_{el\;(x)}}$
(2.6)

Figure 2.3 helps on visualising the different displacements affecting the point mass, starting from a given time ${\textstyle t_{i}}$ until a generic ${\textstyle t}$. The total displacement results from the combination of both motions. Note that, for simplification, only axial elastic deformations are represented.

 Figure 2.3: Representation of the rigid body rotation and axial elastic deformation

Using the rotation matrix obtained in equation 2.1, one can develop the following relations regarding the displacements of the particle:

 ${\displaystyle {\mathbf {U} }\;_{rb\;(x)}\;=\;{\mathbf {{Q}_{(\theta )}} }{\mathbf {r} }_{\,i}\;-\;{\mathbf {r} }_{\,i}=({\mathbf {{Q}_{(\theta )}} }-{\hbox{I}}){\mathbf {r} }_{\,i}\;,\;\;\;\;\;\;\;{\mathbf {U} }\;_{el\;(x)}\;=\;{\mathbf {{Q}_{(\theta )}} }\;{\mathbf {U} }\;_{el\;(\mathbb {\varsigma } )}}$
(2.7)

For the sake of notational simplification, let's simply denote the total displacement in the inertial frame ${\textstyle {\mathbf {U} }\;_{(x)}}$ as ${\textstyle {\mathbf {U} }}$, and the elastic deformations in the corotating frame ${\textstyle {\mathbf {U} }\;_{el\;(\mathbb {\varsigma } )}}$ as ${\textstyle {\mathbf {u} }}$. With this notation, the equations describing the particle's position, velocity and acceleration are:

 ${\displaystyle {\mathbf {U} }={\mathbf {{Q}_{(\theta )}} }(\,{\mathbf {r} }_{\,i}+{\mathbf {u} }\,)-\,{\mathbf {r} }_{\,i}}$
(2.8)
 ${\displaystyle {\mathbf {\dot {U}} }={\mathbf {{Q}_{(\theta )}} }\left[{\widetilde {\Omega }}({\mathbf {r} }_{\,i}+{\mathbf {u} })+{\mathbf {\dot {u}} }\right]}$
(2.9)

 ${\displaystyle {\mathbf {\ddot {U}} }={\mathbf {{Q}_{(\theta )}} }\left[{\underset {Centrifugal}{\underbrace {{\widetilde {\Omega }}^{\,2}({\mathbf {r} }_{\,i}+{\mathbf {u} })} }}+{\underset {Tangential}{\underbrace {{\widetilde {\alpha }}({\mathbf {r} }_{\,i}+{\mathbf {u} })} }}+{\overset {Coriolis}{\overbrace {2{\widetilde {\Omega }}{\mathbf {\dot {u}} }} }}+{\underset {Relative}{\underbrace {\mathbf {\ddot {u}} } }}\right]}$
(2.10)

Where ${\textstyle \alpha }$ stands the angular acceleration set by the driving device. Although until now, the figures accounted only for axial elastic deformations — that is, ${\textstyle {\mathbf {u} }=\left(u\,_{\varsigma },0,0\right)^{T}}$ — the fact is that the previous equations can be used to describe a general deformation in three dimensions: ${\textstyle {\mathbf {u} }=\left(u_{\varsigma },u_{\varrho },u_{\varpi }\right)^{T}}$. From the previous equations, ${\textstyle \Omega }$ and ${\textstyle \alpha }$ are imposed by the driving device, while ${\textstyle {\mathbf {r} }_{\,i}}$ only depends on the initial condition. The elastic kinematic unknowns are, then, ${\textstyle {\mathbf {u} }}$, ${\textstyle {\mathbf {\dot {u}} }}$ and ${\textstyle {\mathbf {\ddot {u}} }}$.

## 2.4 Elasticity notation

Previous to the posing of the equilibrium equations, some definitions regarding the used notation should be clarified. Let ${\textstyle n_{d}}$ denote the number of space dimensions of the problem under consideration, and let ${\textstyle \mho \subset \mathbb {R} ^{n_{d}}}$ be an open set with piecewise smooth boundary ${\textstyle \Gamma }$. A general point, denoted by ${\textstyle {\boldsymbol {x}}}$, is identified with its position vector emanating from the origin. Tensors and vectors will be described in terms of Cartesian components. In indicial notation, differentiation is denoted by:

 ${\displaystyle {\frac {\partial }{\partial x_{i}}}=\partial _{x_{i}}=\partial _{i}\;,\;\;\;\partial _{ii}=\sum _{i=1}^{n_{d}}\partial _{ii}}$
(2.11)

This time, a general form of the elastic problem is to be developed. The displacement field, referred as ${\textstyle {\mathfrak {u}}}$, pretends to represent displacement in a generic case, and should not be confused with the elastic deformations in the corotational frame ${\textstyle {\mathbf {u} }}$, presented in the previous section. Regarding the domain boundary, we shall assume that it does not change it time and admits the following decomposition:

 ${\displaystyle \Gamma =\Gamma _{\mathfrak {u}}^{i}\cup \Gamma _{\sigma }^{i}\,,\;\;\Gamma _{\mathfrak {u}}^{i}\cap \Gamma _{\sigma }^{i}=\varnothing }$
(2.12)

With the previous decomposition, two different boundary conditions can be stated. Regarding the prescribed displacements, the Dirichlet boundary conditions ${\textstyle {\mathfrak {u}}_{i}={\overline {\mathfrak {u}}}^{\,i}}$ are defined on ${\textstyle \Gamma _{\mathfrak {u}}^{i}}$. Note that prescribed displacements may change in time, thus ${\textstyle {\overline {\mathfrak {u}}}_{i}:\Gamma _{u}^{i}\times ]0,T[\rightarrow \mathbb {R} }$, where ${\textstyle ]0,T[}$ is the open time interval of length ${\textstyle T>0}$. On the other hand, regarding the prescribed tractions, the Neumann boundary conditions ${\textstyle {\overline {\mathbf {t} }}^{\,i}=\sigma _{ij}n_{j}}$ are defined on ${\textstyle \Gamma _{\sigma }^{i}}$, where again ${\textstyle {\overline {\mathbf {t} }}_{i}:\Gamma _{\sigma }^{i}\times ]0,T[\rightarrow \mathbb {R} }$.

To complete the definition of the boundary value problem, initial conditions shall be given. Since the governing equation will involve accelerations and thus be second order in time, both initial displacements and velocities must be specified:

 ${\displaystyle {{\mathfrak {u}}_{i}}_{({\boldsymbol {x}},0)}={{\mathfrak {u}}_{i}^{0}}_{({\boldsymbol {x}})}\,,\;\;\;\;{\dot {\mathfrak {u}}}_{i\;({\boldsymbol {x}},0)}={{\dot {\mathfrak {u}}}^{0}}_{i\;({\boldsymbol {x}})}}$
(2.13)

As stated in the previous section, the key is to find the displacements, which are represented by a vector field ${\textstyle {\mathfrak {u}}:\mho \rightarrow \mathbb {R} ^{n_{d}}}$. Using infinitesimal theory (see A.1), the deformation state at a given point ${\textstyle {\boldsymbol {x}}\in \mho }$ is characterised by the infinitesimal strain tensor, which is defined to be the symmetric part of the displacement gradients:

 ${\displaystyle \varepsilon _{ij}:=\left[\nabla ^{s}{\mathfrak {u}}\,\right]_{ij}={\frac {1}{2}}\left({\frac {\partial {u}_{i}}{\partial x_{j}}}+{\frac {\partial {u}_{j}}{\partial x_{i}}}\right)\;\;\;}$
(2.14)

For linear elastostatic systems, the relationship between the Cauchy stress tensor and the strain tensor is given in equation A.6. In the dynamic case, however, the stress at a given time is not only proportional to the strain but also to its derivative (linear viscoelasticity). In indicial form:

 ${\displaystyle \sigma _{ij}=\mathbf {C} _{ijkl}\varepsilon _{kl}+{\overline {\mathbf {D} }}_{ijkl}{\dot {\varepsilon }}_{kl}}$
(2.15)

Where ${\textstyle {\overline {\mathbf {D} }}}$ is the tensor of viscoelasticity coefficients. It is often assumed to be proportional to the elasticity tensor, ${\textstyle {\overline {\mathbf {D} }}={\overline {\beta }}\,\mathbf {C} }$. For the dynamic case, inertia forces — proportional to acceleration — have to be taken into account, and will be treated as body forces per unit volume. Additionally, a damping term proportional to velocity is also to be included:

 ${\displaystyle {\overset {static}{\overbrace {{f}_{i}} }}\rightarrow {\overset {dynamic}{\overbrace {{f}_{i}-\rho {\ddot {\mathfrak {u}}}_{i}-{\overline {\mu }}\;{\dot {\mathfrak {u}}}_{i}} }}}$
(2.16)

Where ${\textstyle \rho :{\overline {\Omega }}\rightarrow \mathbb {R} }$ stands for the material density, and ${\textstyle {\overline {\mu }}}$ is a damping parameter usually defined as ${\textstyle {\overline {\mu }}={\overline {\alpha }}\rho }$, where ${\textstyle {\overline {\alpha }}}$ is known as damping coefficients, often considered to be constant.

## 2.5 Strong form of the elastodynamic equation

Supposing the rotor as perfectly elastic, the equation describing the dynamic behaviour is given by Newton's second law for continuum media, which in indicial notation reads as:

 ${\displaystyle {\frac {\partial \sigma _{ij}}{\partial x_{j}}}+\left({\boldsymbol {f}}_{i}-{\overline {\mu }}\;{\dot {\mathfrak {u}}}_{i}\right)=\rho {\ddot {\mathfrak {u}}}_{i}}$
(2.17)

But to solve the previous equation, the boundary conditions and the elastic model defined in the previous section are required. The strong form of the elastodynamic boundary value problem is then formulated as:

(2.18)

Given ${\textstyle {\boldsymbol {f}}_{i},{\overline {\mathfrak {u}}}^{\,i}\,,\;\;{\overline {\mathbf {t} }}^{\,i}\,,\;\;{\mathfrak {u}}_{i}^{0}}$ and ${\textstyle {\dot {\mathfrak {u}}}_{\,i}^{0}\;}$ , find ${\textstyle {\mathfrak {u}}_{i}:{\overline {\mho }}\times ]0,T[\rightarrow \mathbb {R} \;}$ such that

${\textstyle {\frac {\partial \sigma _{ij}}{\partial x_{j}}}+\left({\boldsymbol {f}}_{i}-{\overline {\mu }}\;{\dot {\mathfrak {u}}}_{i}\right)=\rho \;{\ddot {\mathfrak {u}}}_{i}\;}$ in ${\textstyle \;\mho \times ]0,T[}$

${\textstyle {\mathfrak {u}}_{i}={\overline {\mathfrak {u}}}^{\,i}\;}$ on ${\textstyle \;\Gamma _{\mathfrak {u}}^{i}\times ]0,T[}$

${\textstyle \sigma _{ij}n_{j}={\overline {\mathbf {t} }}^{\,i}\;}$ on ${\textstyle \;\Gamma _{\sigma }^{i}\times ]0,T[}$

${\textstyle {{\mathfrak {u}}_{i}}_{({\boldsymbol {x}},0)}={\mathfrak {u}}_{i\,({\boldsymbol {x}})}^{0}\,}$ and ${\textstyle \;{{\dot {\mathfrak {u}}}_{i\,({\boldsymbol {x}},0)}}={\dot {\mathfrak {u}}}_{i\,({\boldsymbol {x}})}^{0}\;}$ with ${\textstyle \;{\boldsymbol {x}}\in \mho }$

where ${\textstyle \;\sigma _{ij}=\mathbf {C} _{ijkl}\varepsilon _{kl}+{\overline {\mathbf {D} }}_{ijkl}{\dot {\varepsilon }}_{kl}}$

The previous system has no closed form solution, and an analytical solution is only found for simple domains . For more complex geometries and boundary conditions, approximate methods have to be used. There is no need to point out that the effect of rotation does nothing but adds complexity to the equations. In the inertial frame, equation 2.17 yields to the following expression when accounting for the acceleration induced by rotation:

 ${\displaystyle \nabla {\boldsymbol {\sigma }}_{\;{\boldsymbol {x}}}+\;\mathbf {{Q}_{(\theta )}} \left({\mathbf {f} }-{\overline {\mu }}\;{\dot {\mathbf {u} }}\,\right)=\rho \mathbf {{Q}_{(\theta )}} \left[{\widetilde {\Omega }}^{\,2}({\mathbf {r} }_{\,i}+{\mathbf {u} })+{\widetilde {\alpha }}({\mathbf {r} }_{\,i}+{\mathbf {u} })+2{\widetilde {\Omega }}{\mathbf {\dot {u}} }+{\mathbf {\ddot {u}} }\right]}$
(2.19)

### 2.5.1 Transformation of coordinates

One of the problems of the above equation is that the stress field is difficult to pose in terms of the inertial frame, as the elasticity tensor is only known in the corotating coordinates. Using the rotation matrix, stresses can be transformed from one basis to another:

 ${\displaystyle {\sigma }_{\;{\boldsymbol {x}}}=\;\;\mathbf {{Q}_{(\theta )}} \;\;{\sigma }_{\;{\boldsymbol {\varsigma }}}\;\;\mathbf {{Q}_{(\theta )}} ^{\,T}}$
(2.20)

Upon some manipulations, the above expression can be rephrased in Voigt's notation as follows:

 ${\displaystyle \left[{\begin{array}{c}{\sigma _{x}}\\{\sigma _{y}}\\{\sigma _{z}}\\{\tau _{yz}}\\{\tau _{xz}}\\{\tau _{xz}}\end{array}}\right]={\underset {\mathbf {{T}_{(\theta )}} }{\underbrace {\left[{\begin{array}{cccccc}{\cos ^{2}\theta }&{\sin ^{2}\theta }&{0}&{0}&{0}&{-\sin 2\theta }\\{\sin ^{2}\theta }&{\cos ^{2}\theta }&{0}&{0}&{0}&{\sin 2\theta }\\{0}&{0}&{1}&{0}&{0}&{0}\\{0}&{0}&{0}&{\cos \theta }&{\sin \theta }&{0}\\{0}&{0}&{0}&{-\sin \theta }&{\cos \theta }&{0}\\{0.5(\sin 2\theta )}&{-0.5(\sin 2\theta )}&{0}&{0}&{0}&{\cos 2\theta }\end{array}}\right]} }}\left[{\begin{array}{c}{\sigma _{\varsigma }}\\{\sigma _{\varrho }}\\{\sigma _{\varpi }}\\{\tau _{\varrho \varpi }}\\{\tau _{\varsigma \varpi }}\\{\tau _{\varsigma \varrho }}\end{array}}\right]}$
(2.21)

Notice that ${\textstyle \;\mathbf {T} _{(\theta )}^{}\;=\mathbf {T} _{(-\theta )}\;}$.

Then, the relation between the stresses in the inertial ${\textstyle {\boldsymbol {\sigma }}_{\;{\boldsymbol {x}}}}$ and in the corotational frame ${\textstyle {\boldsymbol {\sigma }}_{\;{\boldsymbol {\varsigma }}}}$, and strains in the inertial ${\textstyle {\boldsymbol {\varepsilon }}_{\;{\boldsymbol {x}}}}$ and in the corotational frame ${\textstyle {\boldsymbol {\varepsilon }}_{\;{\boldsymbol {\varsigma }}}}$ in Voigt notation is given by:

 ${\displaystyle {\boldsymbol {\sigma }}_{\;{\boldsymbol {x}}}\;=\;\mathbf {T} _{(\theta )}\;{\boldsymbol {\sigma }}_{\;{\boldsymbol {\varsigma }}}\;,\;\;\;\;\;\;\;{\boldsymbol {\varepsilon }}_{\;{\boldsymbol {x}}}\;=\;\mathbf {T} _{(\theta )}^{T}\;{\boldsymbol {\varepsilon }}_{\;{\boldsymbol {\varsigma }}}}$
(2.22)

An analogous expression can be deduced regarding the elasticity tensor:

 ${\displaystyle {\mathbf {C} }_{\;{\boldsymbol {x}}}\;=\;\mathbf {T} _{(\theta )}\;\;{\mathbf {C} }_{\;{\boldsymbol {\varsigma }}}\;\;\mathbf {T} _{(\theta )}^{T}\;}$
(2.23)

Having defined these relations, equation 2.19 can be rewritten in terms of the known elasticity matrix ${\textstyle {\mathbf {C} }_{\;(\varsigma )}}$:

 ${\displaystyle \nabla \left(\mathbf {T} _{(\theta )}\;{\boldsymbol {\sigma }}_{\;{\boldsymbol {\varsigma }}}\right)+\;\mathbf {{Q}_{(\theta )}} \left({\mathbf {f} }-{\overline {\mu }}\;{\dot {\mathbf {u} }}\,\right)=\rho \mathbf {{Q}_{(\theta )}} \left[{\widetilde {\Omega }}^{\,2}({\mathbf {r} }_{\,i}+{\mathbf {u} })+{\widetilde {\alpha }}({\mathbf {r} }_{\,i}+{\mathbf {u} })+2{\widetilde {\Omega }}{\mathbf {\dot {u}} }+{\mathbf {\ddot {u}} }\right]}$
(2.24)

If the dynamics of the system are studied in the noninertial frame, the above equation reduces to

 ${\displaystyle \nabla {\boldsymbol {\sigma }}_{\;{\boldsymbol {\varsigma }}}+\;\left({\mathbf {f} }-{\overline {\mu }}\;{\dot {\mathbf {u} }}\,\right)=\rho \left[{\widetilde {\Omega }}^{\,2}({\mathbf {r} }_{\,i}+{\mathbf {u} })+{\widetilde {\alpha }}({\mathbf {r} }_{\,i}+{\mathbf {u} })+2{\widetilde {\Omega }}{\mathbf {\dot {u}} }+{\mathbf {\ddot {u}} }\right]}$
(2.25)

### 2.5.2 One-dimensional stationary solution

One of the simplest cases that can be tackled analytically is a one-dimensional stationary approach. In other words, it means that elastic deformations ${\textstyle {\mathbf {u} }}$ are only present in the axial direction ${\textstyle \varsigma }$, while the rotation is considered to be stationary with constant angular velocity ${\textstyle \Omega }$. Moreover, no vibration phenomena is taken into consideration, and so ${\textstyle {\mathbf {\dot {u}} }}$ and ${\textstyle {\mathbf {\ddot {u}} }}$ are neglected. That is to imagine the rotor as a sequence of one-dimensional bars, which can only experience pure traction and compression. For this case where ${\textstyle n_{d}=1}$, neglecting damping and external forces, equation 2.25 yields to

 ${\displaystyle {\frac {dA_{(\varsigma )}\sigma _{\varsigma }}{\partial \varsigma }}=-\rho \,A_{(\varsigma )}\,{\widetilde {\Omega }}^{\,2}({\mathbf {r} }_{\,i}+{\mathbf {u} })\;\;\;\Rightarrow \;\;\;A{\frac {d\left(\mathbf {C} _{\varsigma \varsigma }\;\varepsilon _{\varsigma }\right)}{d\varsigma }}=-\rho A^{\,2}(\varsigma +{{u}_{\varsigma }})}$
(2.26)

Identifying the coefficient ${\textstyle \mathbf {C} _{\varsigma \varsigma }}$ as the Young's Modulus, considered independent on the axial position, and considering constant section ${\textstyle A}$, the above equation can be rewritten as:

 ${\displaystyle E\,A{\frac {d^{2}u_{\varsigma }}{d\varsigma ^{2}}}=-\rho A{\Omega }^{\,2}(\varsigma +{{u}_{\varsigma }})\;\;\;\Rightarrow \;\;\;{\frac {d^{2}u_{\varsigma }}{d\varsigma ^{2}}}=-{\frac {1}{E}}\rho {\Omega }^{\,2}(\varsigma +{{u}_{\varsigma }})}$
(2.27)

Where ${\textstyle \rho }$ stands for the density of the material and ${\textstyle E}$ is the Young's Modulus. The above equation describes the rotation of a 1D structure under constant angular velocity. In the corotational frame, the deformation can be interpreted as the result of applying an axially distributed force, whose modulus is proportional to the radial position (Figure 2.4). Note that the displacement field will not depend on the area of the beam.

The resulting equation can be identified as a second order ordinary differential equation (ODE). Thus, two boundary conditions will be required. Supposing a fixed-free structure, the displacement at the root is to be zero, and the same happens with the stress at the free end, where ${\textstyle \varsigma =R}$ :

 ${\displaystyle \left.{\begin{matrix}u\end{matrix}}\right|_{\varsigma =0}=0\;,\;\;\;\;\left.{\begin{matrix}{\frac {du_{\varsigma }}{d\varsigma }}\end{matrix}}\right|_{\varsigma =R}=0}$
(2.28)

To simplify the solution, let's denote ${\textstyle k={\sqrt {\frac {\rho {\Omega }^{\,2}}{E}}}}$.

 Figure 2.4: 1D static problem interpreted as an axially loaded static bar in the rotating frame

The solutions of the ODE are:

 ${\displaystyle {u_{\varsigma }}_{(\varsigma )}={\frac {\sin(k\varsigma )}{k\cos(kR)}}-\varsigma \;,\;\;\;{\sigma _{\varsigma }}_{(\varsigma )}=E\,\left({\frac {\cos(k\varsigma )}{\cos(kR)}}-1\right)}$
(2.29)

Imagine now an helicopter's rotor, with a constant cross-section and ${\textstyle 4\,m}$ radius rotating at ${\textstyle 60\,rad/s}$. Let's be idealistic, and suppose that the blade can only deform in the axial direction. Moreover, let's imagine that the conditions are perfectly ideal and no source of vibration exists. In that case, the previous equations could be used to compute the deformations and stresses across the beam. For an aluminium beam (${\textstyle E=70\,GPa}$ , ${\textstyle \rho =2700\,kg/m^{3}}$), the obtained results are showed in Figure 2.5, displaying an overall elongation of the blade of ${\textstyle 3\,mm}$:

 Figure 2.5: Analytic solution for an axially loaded bar due to centrifugal forces

For sure, this model cannot be applied to real rotors. However, it is an exact solution to equation 2.25, and will be used latter on this report to validate the developed approximated methods.

# 3. FEM model

Once the strong form of the elastodynamic equation has been posed, is now time to transform it in terms of the finite element method, in such a way that it can be solved using this computation technique. The aim of the present chapter is to pose the weak form of the equation and give a glimpse on how to solve it. It is taken for granted that the reader has basic notions on the topic of FEM, and thus this chapter will not explain the very basics of the method.

## 3.1 Weak form of the elastodynamic equation

For the sake of simplicity, let's take advantage of the formulation presented in the previous chapter while developing the strong form of the problem. To start with, the test and trial functions will be defined. Those are continuous with square-integrable derivative functions, defined over the problem boundary ${\textstyle {\overline {\mho }}}$, with specific boundary values:

 ${\displaystyle {\begin{matrix}\mathbf {\hbox{Test functions }} {\mathcal {}}{V}=\left\{v_{i}:{\overline {\mho }}\rightarrow \mathbb {R} \,,\,\,v_{i}=0{on}\,\Gamma _{u}^{i}\right\}\\\mathbf {\hbox{Trial functions }} {\mathcal {}}{S}=\left\{u_{i}:{\overline {\mho }}\rightarrow \mathbb {R} \,,\,\,u_{i}={\overline {u}}_{i}{on}\,\Gamma _{u}^{i}\right\}\end{matrix}}}$
(3.1)

Denoting as ${\textstyle f_{ij}}$ a ${\textstyle C_{1}}$ function, the Gauss' divergence theorem reads as:

 ${\displaystyle \int _{\mho }{\frac {\partial f_{ij}}{\partial x_{j}}}d\mho =\int _{\Gamma }f_{ij}\;n_{j}d\Gamma }$
(3.2)

The weak form is obtained by integrating over the domain the product of the balance equation (see eq. 2.17) in with the test functions, and by employing the Divergence Theorem to shift the derivatives on the stresses ${\textstyle \sigma _{ij}}$ to the test functions ${\textstyle v_{i}}$. For simplicity, let's denote ${\textstyle (f_{i}-{\overline {\mu }}\;{\dot {\mathfrak {u}}}_{i}-\rho \;{\ddot {\mathfrak {u}}}_{i})}$ as ${\textstyle {f_{i}}}$.

 ${\displaystyle \int _{\mho }v_{i}\left({\frac {\partial \sigma _{ij}}{\partial x_{j}}}+{f_{i}}\right)d\mho \,=\,0\;\;\Rightarrow \;\;\int _{\mho }v_{i}\;{\frac {\partial \sigma _{ij}}{\partial x_{j}}}d\mho +\int _{\mho }v_{i}\;{f_{i}}d\mho \,=\,0}$
(3.3)

Integrating by parts and applying equation 3.2 to the first term:

 ${\displaystyle \int _{\mho }v_{i}{\frac {\partial \sigma _{ij}}{\partial x_{j}}}d\mho =\int _{\Gamma }v_{i}\left(\sigma _{ij}\;n_{j}\right)d\Gamma {-\int }_{\mho }{\frac {\partial v_{i}}{\partial x_{j}}}\;\sigma _{ij}d\mho }$
(3.4)

The first term on the right side of the preceding equation can be expressed as the sum of the integrals over ${\textstyle \Gamma _{u}^{i}{and}\Gamma _{\sigma }^{i}}$. Taking into account that ${\textstyle v_{i}=0{on}\Gamma _{u}^{i}}$, and ${\textstyle \sigma _{ij}n_{j}={\overline {t}}^{\,i}}$ on ${\textstyle \Gamma _{\sigma }^{i}}$ :

 ${\displaystyle \int _{\Gamma }v_{i}\left(\sigma _{ij}\;n_{j}\right)d\Gamma =\sum _{i=1}^{n_{d}}\int _{\Gamma _{u}^{i}}v_{i}\left(\sigma _{ij}\;n_{j}\right)d\Gamma {+\sum }_{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}v_{i}\left(\sigma _{ij}\;n_{j}\right)d\Gamma =\sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}v_{i}\;{\overline {t}}^{\,i}d\Gamma }$
(3.5)

Introducing the latter into equation 3.3:

 ${\displaystyle \sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}v_{i}\;{\overline {t}}^{\,i}d\Gamma {-\int }_{\mho }{\frac {\partial v_{i}}{\partial x_{j}}}\;\sigma _{ij}\;d\mho +\int _{\mho }v_{i}\;{f_{i}}\;d\mho =0}$
(3.6)

As ${\textstyle \sigma _{ij}=\sigma _{ji}}$, the derivatives of the test function can be expressed in terms of the symmetric gradient operator: ${\textstyle {\frac {\partial v_{i}}{\partial x_{j}}}\sigma _{ij}=\left[\nabla ^{s}v\right]_{ij}\sigma _{ij}}$. It will be convenient to express the problem in Voigt's notation (see A.1.2). Thus, the latter expression ${\textstyle \left[\nabla ^{s}v\right]_{ij}\sigma _{ij}}$ transforms into ${\textstyle \left({\boldsymbol {\nabla }}^{s}{\boldsymbol {v}}\right)^{T}{\boldsymbol {\sigma }}}$.

Taking into account the initial conditions defined in the previous chapter, the weak form of the elastodynamic problem using Voigt's notation is written as:

(3.7)

Given ${\textstyle \mathbf {f} ,\;{\overline {\mathfrak {u}}}^{\,i}\,,\;\;{\overline {\mathbf {t} }}^{\,i}\,,\;\;{\mathfrak {u}}^{0}}$ and ${\textstyle {\dot {\mathfrak {u}}}^{0}\;}$ , find ${\textstyle {\mathfrak {u}}_{(t)}\in {\mathcal {}}{S}\;,\;\;t\in [0,T]\;}$ such that

 ${\displaystyle \int _{\mho }\mathbf {v} ^{T}\rho \,{\ddot {\mathfrak {u}}}\;d\mho +\int _{\mho }\mathbf {v} ^{T}\,{\overline {\mu }}\,{\dot {\mathfrak {u}}}\;d\mho +\int _{\mho }\nabla ^{s}\mathbf {v} ^{T}\,{\overline {D}}\,\nabla ^{s}{\dot {\mathfrak {u}}}\,d\mho \,+}$
 ${\displaystyle +\int _{\mho }\nabla ^{s}\mathbf {v} ^{T}\,\mathbf {C} \,\nabla ^{s}{\mathfrak {u}}\,d\mho \;\;=\int _{\mho }{\mathbf {v} }^{T}{\boldsymbol {f}}d\mho +\sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}v_{i}{\overline {t}}^{\,i}d\Gamma \,,\;\forall \mathbf {v} \in {\mathcal {}}{V}}$

${\textstyle {\mathfrak {u}}={\overline {\mathfrak {u}}}^{\,i}\;}$ on ${\textstyle \;\Gamma _{\mathfrak {u}}^{i}\times [0,T]}$

${\textstyle {\mathfrak {u}}_{({\boldsymbol {x}},0)}={\mathfrak {u}}_{({\boldsymbol {x}})}^{0}\,}$ and ${\textstyle \;{\dot {\mathfrak {u}}}_{({\boldsymbol {x}},0)}={\dot {\mathfrak {u}}}_{({\boldsymbol {x}})}^{0}\;}$ with ${\textstyle \;{\boldsymbol {x}}\in \mho }$

The above expression is known as the variational form of the elastodynamic problem, and can be seen as the balance of internal and external virtual works, an equilibrium between stresses and inertia, viscous, body and traction forces. The term ${\textstyle \nabla ^{s}\mathbf {v} ^{T}}$ can be thought as virtual strains, while ${\textstyle \mathbf {v} ^{T}}$ as virtual displacements. The major difference with the strong formulation is that the weaker form of the problem demands less smoothness of the solution, involving only firsts derivatives of the test functions.

The formulation of the weak form of the elastodynamic problem is only the first step for the finite element method to be fully developed. The following sections will introduce the required formulation to tackle the problem in terms of the FEM, allowing its programming.

## 3.2 Global formulation

The key concept of the FEM lies in the domain decomposition into a certain number of elements, finite-sized subdomains. Each element has a given number of nodes, which will depend on the element type. The important fact here is to know that the FEM only computes values of the variable at the nodes. Then, the values of non-nodal points are approximated by interpolation of the nodal values. The FEM equations are formulated in such a way that continuity of the field variables is ensured. However, interelement continuity of the gradients of the variable does not often exist. The latter is a problem affecting some variables of interest such as strains and stresses, and is solved by increasing the number of elements used.

Following with the concept of domain decomposition, it is said that the domain ${\textstyle {\overline {\mho }}}$ is discretised into element domains ${\textstyle {\overline {\mho }}\,^{e}}$. The geometric shape resulting from the domain division is known as mesh. If the domain can be thought as a polygon, and ${\textstyle n_{el}}$ is the total number of elements in which the domain is split, the decomposition can be written as

 ${\displaystyle {\overline {\mho }}={\overline {\mho }}^{\,1}\cup {\overline {\mho }}^{\,2}\cdots \cup {\overline {\mho }}^{\,n_{el}}=\bigcup _{e=1}^{n_{el}}{\overline {\mho }}^{\,e}}$
(3.8)

If ${\textstyle n_{pt}}$ is the total number of nodes of the domain, then the FEM has to compute the field variable (in this case, displacement) in each of those nodes. For ${\textstyle n_{d}}$ dimensions, each nodal displacement will be a ${\textstyle n_{d}}$-sized vector. We can define a global vector of displacements d${\textstyle \,\in \,\mathbb {R} ^{n_{d}n_{pt}}}$ and variations c${\textstyle \,\in \,\mathbb {R} ^{n_{d}n_{pt}}}$ :

 ${\displaystyle \mathbf {d} :=\left[{\begin{array}{c}{\mathbf {d} _{1}}\\{\mathbf {d} _{2}}\\{\vdots }\\{\mathbf {d} _{\,n_{pt}}}\end{array}}\right],\quad \mathbf {c} :=\left[{\begin{array}{c}{\mathbf {c} _{1}}\\{\mathbf {c} _{2}}\\{\vdots }\\{\mathbf {c} _{\,n_{pt}}}\end{array}}\right]}$
(3.9)

It is important to notice that not all the components of the previous vectors are unknown: ${\textstyle \mathbf {r} \subset \left\{1,2\ldots n_{d}n_{pt}\right\}}$ is defined as the set of constrained degrees of freedom (DOF), along which displacement is known. Thus, ${\textstyle \mathbf {d} _{r}={\overline {\mathfrak {u}}}}$ and ${\textstyle \mathbf {c} _{r}=0}$. On the other hand, ${\textstyle \mathbf {l} }$ is known as the set of unconstrained degrees of freedom, the global DOFs along which displacement is unknown.

As it has been stated in the first paragraph of this section, the solution along non-nodal points is interpolated using the nodal values. The interpolation functions are known as shape functions, and in the Galerkin method the same functions are used to interpolate both shape and trial functions. Let's define the continuous and with piecewise continuous derivative scalar shape function ${\textstyle N_{A}:{\overline {\mho }}\rightarrow \mathbb {R} }$, where ${\textstyle 1\leq A\leq n_{pt}}$.

In the FEM, nodal shape functions ${\textstyle N_{A}}$ are defined in such a way that their value is ${\textstyle \mathbf {1} }$ at node ${\textstyle A}$, and ${\textstyle \mathbf {0} }$ at every other node. In problems tackling with vector fields, a matrix of nodal shape functions is defined. Using the same function to interpolate displacements in all dimensions, the nodal matrix of shape functions for a three-dimensional case is:

 ${\displaystyle \mathbf {N} _{A}:=N_{A}\,\mathbf {I} =N_{A}\left[{\begin{array}{ccc}{1}&{0}&{0}\\{0}&{1}&{0}\\{0}&{0}&{1}\end{array}}\right]\,,\;\;\mathbf {N} _{A}\in \mathbb {R} ^{n_{d}\times n_{d}}}$
(3.10)

A global matrix of shape functions N ${\textstyle \rightarrow \mathbb {R} ^{n_{d}\times n_{d}\,n_{pt}}}$ can be defined, involving all the nodal matrices: ${\textstyle \mathbf {N} :=\left[{\begin{array}{llll}{\mathbf {N} _{1}}&{\mathbf {N} _{2}}&{\cdots }&{\mathbf {N} _{n_{pt}}}\end{array}}\right]}$. In the same line, a global matrix ${\textstyle \mathbf {B} }$ involving the gradient of the shape functions is defined: ${\textstyle \mathbf {B} ={\boldsymbol {\nabla }}^{s}\,\mathbf {N} }$.

The next step is to construct finite-dimensional approximations of the test and trial functions, denoted as ${\textstyle {\mathcal {}}{V}^{h}}$ and ${\textstyle {\mathcal {}}{S}^{h}}$, respectively. If one thinks of these approximations as subsets of ${\textstyle {\mathcal {}}{V}}$ and ${\textstyle {\mathcal {}}{S}}$, and if ${\textstyle u^{h}\subset {\mathcal {}}{S}^{h}}$ and ${\textstyle v^{h}\subset {\mathcal {}}{V}^{h}}$, then the approximated functions will meet the boundary conditions defined in 3.1. The finite dimensional space of test and trial functions is defined as the space including all the functions of the form

 ${\displaystyle {\begin{matrix}{\mathcal {}}{V}^{h}\rightarrow \mathbf {v} _{\,({\boldsymbol {x}})}^{h}=\mathbf {N} _{\,({\boldsymbol {x}})}\,\mathbf {c} ,{with}\mathbf {c} _{\mathrm {r} }=0\\\\{\mathcal {}}{S}^{h}\rightarrow {\mathfrak {u}}_{\,({\boldsymbol {x}},\,{t})}^{h}=\mathbf {N} _{\,({\boldsymbol {x}})}\,\mathbf {d} _{(t)},{with}{\mathbf {d} _{\mathrm {r} }}_{(t)}={\overline {\mathfrak {u}}}_{(t)}\end{matrix}}}$
(3.11)

In a similar way, velocity ${\textstyle {\dot {\mathfrak {u}}}}$ and acceleration ${\textstyle {\ddot {\mathfrak {u}}}}$ are approximated as ${\textstyle {\dot {\mathfrak {u}}}^{h}=\mathbf {N} _{\,({\boldsymbol {x}})}\,{\dot {\mathbf {d} }}_{(t)}}$ and ${\textstyle {\ddot {\mathfrak {u}}}^{h}=\mathbf {N} _{\,({\boldsymbol {x}})}\,{\ddot {\mathbf {d} }}_{(t)}}$, respectively. With all the previous definitions being formulated, the approximated weak form of the problem can be posed:

(3.12)

Given ${\textstyle \mathbf {f} ,\;{\overline {\mathfrak {u}}}^{\,i}\,,\;\;{\overline {\mathbf {t} }}^{\,i}\,,\;\;{\mathfrak {u}}^{0}}$ and ${\textstyle {\dot {\mathfrak {u}}}^{0}\;}$ , find ${\textstyle {\mathfrak {u}}_{(t)}^{h}\in {\mathcal {}}{S}^{h}\;,\;\;t\in [0,T]\;}$ such that

 ${\displaystyle \int _{\mho }{\mathbf {v} \,^{h}}^{T}\rho \,{\ddot {\mathfrak {u}}}^{h}d\mho +\int _{\mho }{\mathbf {v} \,^{h}}^{T}\,{\overline {\mu }}\,{\dot {\mathfrak {u}}}^{h}\,d\mho +\int _{\mho }\nabla ^{s}{\mathbf {v} \,^{h}}^{T}\,{\overline {D}}\,\nabla ^{s}{\dot {\mathfrak {u}}}^{h}\,d\mho \,+}$
 ${\displaystyle +\int _{\mho }\nabla ^{s}{\mathbf {v} \,^{h}}^{T}\,\mathbf {C} \,\nabla ^{s}{\mathfrak {u}}^{h}\,d\mho \;\;=\int _{\mho }{\mathbf {v} \,^{h}}^{T}{\boldsymbol {f}}d\mho +\sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}v_{i}^{h}{\overline {t}}^{\,i}d\Gamma \,,\;\forall \mathbf {v} \in {\mathcal {}}{V}}$

${\textstyle {\mathfrak {u}}^{h}={\overline {\mathfrak {u}}}^{\,i}\;}$ on ${\textstyle \;\Gamma _{\mathfrak {u}}^{i}\times [0,T]}$

${\textstyle {\mathfrak {u}}_{({\boldsymbol {x}},0)}^{h}={\mathfrak {u}}_{({\boldsymbol {x}})}^{0}\,}$ and ${\textstyle \;{\dot {\mathfrak {u}}}_{({\boldsymbol {x}},0)}^{h}={{\dot {\mathfrak {u}}}^{0}}_{({\boldsymbol {x}})}\;}$ with ${\textstyle \;{\boldsymbol {x}}\in \mho }$

### 3.2.1 Approximated rotating elastodynamic model

Until now, rotation has not been taken into account when modelling the weak form of the elastodynamic problem. Introducing the kinematic equations obtained in section 2.3, substituting the approximations presented in equation 3.11 and using the algebraic property ${\textstyle (\mathbf {AB} )^{\mathrm {T} }=\mathbf {B} ^{\mathrm {T} }\mathbf {A} ^{\mathrm {T} }}$, the above expression can be written in the noninertial frame as

 ${\displaystyle \mathbf {c} ^{T}\left\{\;{\overset {\mathbf {M} }{\overbrace {\left(\int _{\mho }^{}\mathbf {N} ^{T}\rho \;\mathbf {N} \;d\mho \right)} }}\;{\ddot {\mathbf {d} }}\;+\;{\overset {\mathbf {D} }{\overbrace {\left({\underset {\mathbf {D} _{\,static}}{\underbrace {\int _{\mho }^{}\mathbf {N} ^{T}{\overline {\mu }}\;\mathbf {N} \;d\mho \;+\;\int _{\mho }^{}\mathbf {B} ^{T}{\overline {\mathbf {D} }}\;\mathbf {B} \;d\mho } }}\;+\;{\underset {\mathbf {D} _{\,rot}}{\underbrace {2\int _{\mho }^{}\mathbf {N} ^{T}\rho \;{\widetilde {\Omega }}\;\mathbf {N} \;d\mho } }}\right)} }}\;{\dot {\mathbf {d} }}\;+\right.}$
 ${\displaystyle +\;{\overset {\mathbf {K} }{\overbrace {\left(\;{\underset {\mathbf {K} _{\,static}}{\underbrace {\int _{\mho }^{}\mathbf {B} ^{T}\mathbf {C} \;\mathbf {B} \;d\mho } }}\;+\;{\underset {\mathbf {K} _{\,rot}}{\underbrace {\int _{\mho }^{}\mathbf {N} ^{T}\rho \;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\;\mathbf {N} \;d\mho } }}\right)} }}\;{\mathbf {d} }\;-}$

 ${\displaystyle \left.\;-\;{\overset {\mathbf {F} }{\overbrace {\left(\;{\underset {\mathbf {F} _{\,static}}{\underbrace {\int _{\mho }{\mathbf {N} }^{T}{\boldsymbol {f}}d\mho +\sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}{\overline {\mathbf {N} }}_{i}^{T}\;{\overline {t}}^{\,i}d\Gamma } }}\;-\;{\underset {\mathbf {F} _{\,rot}}{\underbrace {\int _{\mho }^{}\mathbf {N} ^{T}\rho \;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\;\mathbf {r} _{i}\;d\mho } }}\right)} }}\right\}=0}$
(3.13)

In compact form:

 ${\displaystyle \mathbf {c} ^{T}\left[\mathbf {M} \,{\ddot {\mathbf {d} }}\;+\;\mathbf {D} \,{\dot {\mathbf {d} }}\;+\;\mathbf {K} \,{\mathbf {d} }-\mathbf {F} \right]=0}$

 ${\displaystyle \mathbf {c} ^{T}\left[\mathbf {M} \,{\ddot {\mathbf {d} }}\;+\;(\mathbf {D} _{\;static}+\mathbf {D} _{\,rot})\,{\dot {\mathbf {d} }}\;+\;(\mathbf {K} _{\,static}+\mathbf {K} _{\,rot})\,{\mathbf {d} }-(\mathbf {F} _{\,static}-\mathbf {F} _{\,rot})\right]=0}$
(3.14)

Where the subscript "static" refers to those matrices that are general in the FEM when applied to static structures, whereas those matrices with the subscript "rot" are particular of this specific case of rotating structure. Recall that in the above equation, d stands for the approximate elastic displacement in the rotative frame — ${\textstyle \,\mathbf {u} _{(\,{\boldsymbol {x}},t)}^{h}=\mathbf {N} _{(x)}\;\mathbf {d} _{(t)}\,}$.

The obtained matrices are very common in physical systems. ${\textstyle \mathbf {M} }$ is known as mass matrix, ${\textstyle \mathbf {D} }$ as damping matrix, ${\textstyle \mathbf {K} }$ as stiffness matrix and ${\textstyle \mathbf {F} }$ as force vector. Quite often, ${\textstyle {\overline {\mathbf {D} }}}$ and ${\textstyle {\overline {\boldsymbol {\mu }}}}$ are not known. A common approach is to suppose that the system is governed by the Rayleigh damping, which assumes ${\textstyle {\overline {\boldsymbol {\mu }}}={\overline {\alpha }}\rho }$ and ${\textstyle {\overline {\mathbf {D} }}={\overline {\beta }}\mathbf {C} }$, where ${\textstyle {\overline {\alpha }}}$ and ${\textstyle {\overline {\beta }}}$ are assumed to be constant parameters. Thus, the static damping matrix in the noninertial frame can be expressed as:

 ${\displaystyle \mathbf {D} \;_{static}={\overline {\alpha }}\,\mathbf {M} +{\overline {\beta }}\,\mathbf {K} \;_{static}}$
(3.15)

It is important to notice that, in the case where the angular velocity ${\textstyle \Omega }$ is constant, the stiffness and damping matrices and the force vector due to rotation can be written as:

 ${\displaystyle \mathbf {K} _{\,rot}=\int _{\mho }^{}\mathbf {N} ^{T}\rho \;{\widetilde {\Omega }}^{2}\;\mathbf {N} \;d\mho =-{\Omega }^{2}\int _{\mho }^{}\mathbf {N} ^{T}\rho \;{\begin{bmatrix}1&0&0\\0&1&0\\0&0&0\end{bmatrix}}\;\mathbf {N} \;d\mho }$
 ${\displaystyle \mathbf {D} _{\,rot}=2\int _{\mho }^{}\mathbf {N} ^{T}\rho \;{\widetilde {\Omega }}\;\mathbf {N} \;d\mho =2{\Omega }\int _{\mho }^{}\mathbf {N} ^{T}\rho \;{\begin{bmatrix}0&-1&0\\1&0&0\\0&0&0\end{bmatrix}}\;\mathbf {N} \;d\mho }$
 ${\displaystyle \mathbf {F} _{\,rot}=\int _{\mho }^{}\mathbf {N} ^{T}\rho \;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\;\mathbf {r} _{i}\;d\mho =-{\Omega }^{2}\int _{\mho }^{}\mathbf {N} ^{T}\rho \;{\begin{bmatrix}1&0&0\\0&1&0\\0&0&0\end{bmatrix}}\;\mathbf {r} _{i}\;d\mho }$
(3.16)

In this case, all the above expressions are constant in time. The key fact here is that, for the two-dimensional case, ${\textstyle \mathbf {K} _{\,rot}}$ is proportional to the mass matrix by a factor of ${\textstyle -\Omega ^{2}}$. Then:

 ${\displaystyle \mathbf {K} =\mathbf {K} _{\,static}-\Omega ^{2}\mathbf {M} }$
(3.17)

This is a softening effect induced by the non-inertial character of the frame reference where the equations are posed. This behaviour has already been reviewed in the introductory part of this report (see section 1.2.6). This softening causes a reduction in the natural frequencies of the structure, compared to those of the static case. However, this effect is rarely seen in actual rotors. In reality, softening is a phenomena eclipsed by the so-called stiffening effect. The problem is that simple rotor models do not account for the latter, and thus the overall effect is of softening [3].

For sure, equation 3.14 can be posed in the inertial frame of reference using the rotation matrix and applying the transformation of stresses from equation 2.23:

 ${\displaystyle \mathbf {M} _{\;inertial}=\int _{\mho }^{}\mathbf {N} ^{T}\mathbf {{Q}_{(\theta )}} \;\rho \;\mathbf {N} \;d\mho }$
 ${\displaystyle \mathbf {D} _{\;inertial}=\int _{\mho }^{}\mathbf {N} ^{T}\mathbf {{Q}_{(\theta )}} \;{\overline {\mu }}\;\mathbf {N} \;d\mho \;+\;\int _{\mho }^{}\mathbf {B} ^{T}\mathbf {T} _{(\theta )}\;{\overline {\mathbf {D} }}\;\mathbf {T} _{(\theta )}^{T}\;\mathbf {B} \;d\mho \;+\;2\int _{\mho }^{}\;\mathbf {N} ^{T}\mathbf {{Q}_{(\theta )}} \;\rho \;{\widetilde {\Omega }}\;\mathbf {N} \;d\mho }$
 ${\displaystyle \mathbf {K} _{\;inertial}=\int _{\mho }^{}\mathbf {B} ^{T}\mathbf {T} _{(\theta )}\;\mathbf {C} \;\mathbf {T} _{(\theta )}^{T}\;\mathbf {B} \;d\mho \;+\;\int _{\mho }^{}\mathbf {N} ^{T}\mathbf {{Q}_{(\theta )}} \;\rho \;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\;\mathbf {N} \;d\mho }$
 ${\displaystyle \mathbf {F} _{\;inertial}=\int _{\mho }{\mathbf {N} }^{T}\mathbf {{Q}_{(\theta )}} \;{\boldsymbol {f}}d\mho +\sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}\mathbf {N} _{i}^{T}\;\mathbf {T} _{(\theta )}\;{\overline {t}}^{\,i}d\Gamma \;-\;{\int _{\mho }^{}\mathbf {N} ^{T}\mathbf {{Q}_{(\theta )}} \;\rho \;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\;\mathbf {r} _{i}\;d\mho }}$
(3.18)

### 3.2.2 Block matrix decomposition

The above equations cannot be solved unless boundary values of displacement are specified. Moreover, the involved matrices do not admit inverse, so the only way to tackle the system is to break the presented matrices into block matrices. Remember that ${\textstyle \mathbf {r} }$ stands for the set of DOFs where displacement is prescribed, whereas ${\textstyle \mathbf {l} }$ holds for the set of DOFs where the displacement field is to be computed. Hence, naming ${\textstyle \mathbf {J} _{rl}}$ the matrix made-up by the elements forming the intersection of the ${\textstyle \mathbf {r} }$ rows and ${\textstyle \mathbf {l} }$ columns of ${\textstyle \mathbf {J} }$, equation 3.14 is split into:

 ${\displaystyle {\begin{bmatrix}{\mathbf {c} _{\mathbf {l} }}^{T}&{\mathbf {c} _{\mathbf {r} }}^{T}\end{bmatrix}}\left({\begin{bmatrix}\mathbf {M} _{\mathbf {ll} }&\mathbf {M} _{\mathbf {lr} }\\\mathbf {M} _{\mathbf {rl} }&\mathbf {M} _{\mathbf {rr} }\end{bmatrix}}{\begin{bmatrix}{\ddot {\mathbf {d} }}_{\mathbf {l} }\\{\ddot {\mathbf {d} }}_{\mathbf {r} }\end{bmatrix}}+{\begin{bmatrix}\mathbf {D} _{\mathbf {ll} }&\mathbf {D} _{\mathbf {lr} }\\\mathbf {D} _{\mathbf {rl} }&\mathbf {D} _{\mathbf {rr} }\end{bmatrix}}{\begin{bmatrix}{\dot {\mathbf {d} }}_{\mathbf {l} }\\{\dot {\mathbf {d} }}_{\mathbf {r} }\end{bmatrix}}+{\begin{bmatrix}\mathbf {K} _{\mathbf {ll} }&\mathbf {K} _{\mathbf {lr} }\\\mathbf {K} _{\mathbf {rl} }&\mathbf {K} _{\mathbf {rr} }\end{bmatrix}}{\begin{bmatrix}{\mathbf {d} }_{\mathbf {l} }\\{\mathbf {d} }_{\mathbf {r} }\end{bmatrix}}-{\begin{bmatrix}\mathbf {F} _{\mathbf {l} }\\\mathbf {F} _{\mathbf {r} }+\mathbf {R} \end{bmatrix}}\right)=0}$
(3.19)

Where ${\textstyle \mathbf {R} }$ stands for the elastic reactions at the nodes with prescribed displacement, which are unknown. The above expression can be thought as a system of two equations:

 ${\displaystyle {\begin{matrix}\mathbf {c} _{\mathbf {l} }^{T}\left(\mathbf {M} _{\mathbf {ll} }{\ddot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {M} _{\mathbf {lr} }{\ddot {\mathbf {d} }}_{\mathbf {r} }\;+\;\mathbf {D} _{\mathbf {ll} }{\dot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {D} _{\mathbf {lr} }{\dot {\mathbf {d} }}_{\mathbf {r} }\;+\;\mathbf {K} _{\mathbf {ll} }{\mathbf {d} }_{\mathbf {l} }\;+\;\mathbf {K} _{\mathbf {lr} }{\mathbf {d} }_{\mathbf {r} }\;-\;\mathbf {F} _{\mathbf {l} }\right)=0\\\\\mathbf {c} _{\mathbf {r} }^{T}\left(\mathbf {M} _{\mathbf {rl} }{\ddot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {M} _{\mathbf {rr} }{\ddot {\mathbf {d} }}_{\mathbf {r} }\;+\;\mathbf {D} _{\mathbf {rl} }{\dot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {D} _{\mathbf {rr} }{\dot {\mathbf {d} }}_{\mathbf {r} }\;+\;\mathbf {K} _{\mathbf {rl} }{\mathbf {d} }_{\mathbf {l} }\;+\;\mathbf {K} _{\mathbf {rr} }{\mathbf {d} }_{\mathbf {r} }\;-\;\mathbf {F} _{\mathbf {r} }\;-\;\mathbf {R} \right)=0\end{matrix}}}$
(3.20)

The above equations are fulfilled for all ${\textstyle \mathbf {c} ^{T}}$ if and only if the terms multiplying both ${\textstyle \mathbf {c} _{\mathbf {l} }^{T}}$ and ${\textstyle \mathbf {c} _{\mathbf {r} }^{T}}$ vanish. Solving for the unknown displacements and reactions, and substituting ${\textstyle {\ddot {\mathbf {d} }}_{\mathbf {r} }}$, ${\textstyle {\dot {\mathbf {d} }}_{\mathbf {r} }}$ and ${\textstyle {\mathbf {d} }_{\mathbf {r} }}$ by the corresponding boundary values ${\textstyle {\overline {\ddot {\mathbf {u} }}}}$, ${\textstyle {\overline {\dot {\mathbf {u} }}}}$ and ${\textstyle {\overline {\mathbf {u} }}}$, respectively, the resulting equations are:

 ${\displaystyle {\begin{matrix}\mathbf {M} _{\mathbf {ll} }{\ddot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {D} _{\mathbf {ll} }{\dot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {K} _{\mathbf {ll} }{\mathbf {d} }_{\mathbf {l} }=\mathbf {F} _{\mathbf {l} }-(\;\mathbf {M} _{\mathbf {lr} }{\overline {\ddot {\mathbf {u} }}}\;+\;\mathbf {D} _{\mathbf {lr} }{\overline {\dot {\mathbf {u} }}}\;+\;\mathbf {K} _{\mathbf {lr} }{\overline {\mathbf {u} }})\\\\\mathbf {R} =\mathbf {M} _{\mathbf {rl} }{\ddot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {D} _{\mathbf {rl} }{\dot {\mathbf {d} }}_{\mathbf {l} }\;+\;\mathbf {K} _{\mathbf {rl} }{\mathbf {d} }_{\mathbf {l} }\;+\;\mathbf {M} _{\mathbf {rr} }{\overline {\ddot {\mathbf {u} }}}\;+\;\mathbf {D} _{\mathbf {rr} }{\overline {\dot {\mathbf {u} }}}\;+\;\mathbf {K} _{\mathbf {rr} }{\overline {\mathbf {u} }}\;-\;\mathbf {F} _{\mathbf {r} }\end{matrix}}}$
(3.21)

### 3.2.3 Reaction computation

Notice that the reactions can only be found once the displacements in the unconstrained degrees of freedom are known. If the equations are posed in the noninertial frame of reference, the obtained value of reaction is the expected from a static structure under the effect of a certain deformation. To account for the total reaction, the rigid body moment associated to the rotation motion is to be known:

 ${\displaystyle \mathbf {W} _{rb}=\mathbf {J} _{\mathbf {z} }\;\alpha }$
(3.22)

Where ${\textstyle \mathbf {W} _{rb}}$ stands for the reaction torque associated to the rigid body motion, ${\textstyle \mathbf {J} _{\mathbf {z} }}$ is the inertia moment around the rotation axis and ${\textstyle \alpha }$ is the angular acceleration. If this reaction torque is added to the torque generated by the elastic reactions in the noninertial frame ${\textstyle \mathbf {R} }$, the overall torque is obtained. In order to find the value of ${\textstyle \mathbf {W} _{rb}}$, it will be necessary to compute ${\textstyle \mathbf {J} _{\mathbf {z} }}$, which fortunately can be posed in terms of the FEM. Defining the rigid body modes ${\textstyle {\mathcal {}}{R}}$ in three dimensions:

 ${\displaystyle {\mathcal {}}{R}={\begin{bmatrix}1&0&0&0&\Delta z_{1}&-\Delta y_{1}\\0&1&0&-\Delta z_{1}&0&\Delta x_{1}\\0&0&1&\Delta y_{1}&-\Delta x_{1}&0\\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots \\1&0&0&0&\Delta z_{n_{pt}}&-\Delta y_{n_{pt}}\\0&1&0&-\Delta z_{n_{pt}}&0&\Delta x_{n_{pt}}\\0&0&1&\Delta y_{n_{pt}}&-\Delta x_{n_{pt}}&0\\\end{bmatrix}}={\begin{bmatrix}\mathbf {I} &-\Delta {\widetilde {\mathbf {r} }}\\\vdots &\vdots \end{bmatrix}}}$
(3.23)

Where ${\textstyle {\mathcal {}}{R}}$ is a matrix with 6 columns, corresponding to the rigid body degrees of freedom: three translations and three rotations. In the above expressions, the first three columns are the translation modes, and can be represented by the identity matrix I. On the other hand, the last three columns are the rotation modes, and are computed as the spin of the distance with respect to the rotation axis ${\textstyle -\Delta {\widetilde {\mathbf {r} }}}$. The rigid body modes matrix has as many rows as degrees of freedom, which means that the sequence presented above is repeated for each node. In two-dimensional problems, only one rotational and two translational modes exist.

The key fact is that ${\textstyle {\mathcal {}}{R}}$ can be used to normalise the mass matrix ${\textstyle \mathbf {M} }$ in such a way that information about the structure is obtained. For a generic two-dimensional case it is found that:

 ${\displaystyle {\underset {3\times n_{pt}}{\underbrace {{\mathcal {}}{R}^{T}} }}\;{\underset {n_{pt}\times n_{pt}}{\underbrace {\mathbf {M} } }}\;{\underset {n_{pt}\times 3}{\underbrace {{\mathcal {}}{R}} }}={\begin{bmatrix}\mathbf {m} &0&0\\0&\mathbf {m} &0\\0&0&\mathbf {J} _{\mathbf {z} }\end{bmatrix}}}$
(3.24)

Where m stands for the total mass of the structure, and ${\textstyle \mathbf {J} _{\mathbf {z} }}$ is the inertia moment around the rotation axis. These properties may seem to be trivial to compute, but determining either the mass or the inertia can be a difficult task if the structure is not defined by a simple geometry. Moreover, when the geometry is simple and those values can be computed analytically, they can be used to asses whether ${\textstyle \mathbf {M} }$ has been computed correctly by comparing the analytic values of mass and inertia with those obtained using the above numerical expression.

Once the inertia is known, the torque ${\textstyle \mathbf {W} _{rb}}$ associated to the rigid body motion can be computed. Then, the reaction torque due to elastic reactions ${\textstyle \mathbf {W} _{el}}$ is to be computed by multiplying each term of ${\textstyle \mathbf {R} }$ by the distance with respect to the rotation axis. For a generic case, the elastic reaction torque is computed in terms of each of the components of ${\textstyle \mathbf {R} }$:

 ${\displaystyle \mathbf {w} _{el\;}={\begin{bmatrix}0&-\Delta z\,_{1}&\Delta y\,_{1}\\\Delta z\,_{1}&0&-\Delta x\,_{1}\\-\Delta y\,_{1}&\Delta x\,_{1}&0\\\vdots &\vdots &\vdots \\0&-\Delta z_{{n}_{pt}}&\Delta y_{{n}_{pt}}\\\Delta z_{{n}_{pt}}&0&-\Delta x_{{n}_{pt}}\\-\Delta y_{{n}_{pt}}&\Delta x_{{n}_{pt}}&0\\\end{bmatrix}}{\begin{pmatrix}{R_{x}}\,_{1}\\{R_{y}}\,_{2}\\{R_{z}}\,_{3}\\\vdots \\{R_{x}}_{{n}_{pt}}\\{R_{y}}_{{n}_{pt}}\\{R_{z}}_{{n}_{pt}}\end{pmatrix}}={\begin{bmatrix}\Delta {\widetilde {\mathbf {r} }}\\\vdots \end{bmatrix}}\mathbf {R} }$
(3.25)

The total reaction torque ${\textstyle \mathbf {W} _{el\;}}$ is an ${\textstyle n_{d}}$ vector computed by adding the contributions of ${\textstyle \mathbf {w} _{el\;}}$ for every spatial direction. Finally, the total torque is:

 ${\displaystyle \mathbf {T} =\mathbf {W} _{rb}+\mathbf {W} _{el}}$
(3.26)

When the kinematics of the problem are known — which is the case — one of the the most important unknowns to determine are the total forces and moments acting on the structure. The value of ${\textstyle \mathbf {T} }$ gives a hint on how much work needs to be injected into the blade to ensure rotation at the given angular speed and acceleration.

## 3.3 Local formulation

The integrals presented in the previous sections cannot be computed at once unless the geometry of the domain is tremendously simple, which is not the case. The grace of the FEM is that breaks the domain into smaller elements with well defined geometric properties, where the previous integrals are approximated using Gauss Quadrature, allowing computer implementation . In order to achieve this goal, a new local frame within a single element ${\textstyle \mho ^{e}}$ is defined. Let's assume the model consists of ${\textstyle n_{el}}$ elements, and take ${\textstyle e}$ as the variable index for the elements ${\textstyle 1\leq e\leq n_{el}}$.

Let's introduce ${\textstyle {\boldsymbol {x}}^{e}\;\in {\overline {\mho }}^{\boldsymbol {\,e}}}$ , the vector of global or physical coordinates in the local frame, containing the coordinates of a given element ${\textstyle \mathbf {e} }$. Following this local approach, let's define ${\textstyle {\boldsymbol {\xi }}=(\xi \;\eta \;\zeta \;)^{T}\;\in {\overline {\mho }}_{\boldsymbol {\xi }}\;}$ as the vector of local or element coordinates in the parent domain. In order to visualise the mapping from the parent domain to the physical domain using isoparametric linear elements, let's consider a simple 1D case. The transformation from one domain to another is represented in Figure 3.1.

 Figure 3.1: Representation of the mapping from the physical to the parent domain for a linear 1D element

The same concept can be applied to 2D and 3D elements. For a generic case, the mapping is constructed using the same shape functions employed in the interpolation of ${\textstyle {\mathfrak {u}}^{h}}$ (isoparametric elements). Defining ${\textstyle n^{e}}$ as the number of nodes of the element ${\textstyle e}$, the mapping is expressed as:

 ${\displaystyle {\boldsymbol {x}}_{(\xi )}^{e}={\begin{bmatrix}x_{1}^{e}&x_{2}^{e}&\cdots &x_{n_{e}}^{e}\\y_{1}^{e}&y_{2}^{e}&\cdots &y_{n_{e}}^{e}\\z_{1}^{e}&z_{2}^{e}&\cdots &z_{n_{e}}^{e}\end{bmatrix}}{\begin{bmatrix}{N}_{1}^{e}\\{N}_{2}^{e}\\\vdots \\{N}_{n_{e}}^{e}\end{bmatrix}}=\mathbf {X} ^{e}{{\mathcal {}}{N}_{(\xi )}^{e}}^{T}}$
(3.27)

Were ${\textstyle \mathbf {X} ^{e}}$ stands for the nodal coordinates in the physical domain, and ${\textstyle {\mathcal {}}{N}_{(\xi )}^{e}\in \mathbb {R} ^{1\times n^{e}}}$ is the matrix of shape functions for scalar valued fields. In a similar way, one may define the matrix of shape functions for vector valued fields ${\textstyle \mathbf {N} ^{e}=[{N}_{1}^{e}\;\mathbf {I} \;\;\;{N}_{2}^{e}\;\mathbf {I} \;\;\;\cdots \;\;\;{N}_{n_{e}}^{e}\;\mathbf {I} ]\in \mathbb {R} ^{n_{d}\times n_{d}\;n^{e}}}$ , where ${\textstyle \mathbf {I} }$ is the ${\textstyle n_{d}\times n_{d}}$ identity matrix.

The mapping can be interpreted as a change of coordinates, and thus it is a transformation with an associated Jacobian matrix ${\textstyle \mathbf {J} ^{e}}$. It allows the differential mapping between physical and parent domain, and vice versa:

 $\begin{array}{c}\left[\begin{array}{l}dx\\ dy\\ dz\end{array}\right]=\left[\begin{array}{lll}\frac{\mathrm{\partial }x}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }x}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }x}{\mathrm{\partial }\zeta }\\ \frac{\mathrm{\partial }y}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }y}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }y}{\mathrm{\partial }\zeta }\\ \frac{\mathrm{\partial }z}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }z}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }z}{\mathrm{\partial }\zeta }\end{array}\right]\left[\begin{array}{l}d\xi \\ d\eta \\ d\zeta \end{array}\right]\phantom{\rule{1em}{0ex}}⇒\phantom{\rule{1em}{0ex}}d\mathbit{x}={\mathbf{J}}^{e}d\mathbit{\xi }\phantom{\rule{thickmathspace}{0ex}}\phantom{\rule{1em}{0ex}}\\ \\ \left[\begin{array}{l}d\xi \\ d\eta \\ d\zeta \end{array}\right]=\left[\begin{array}{lll}\frac{\mathrm{\partial }\xi }{\mathrm{\partial }x}& \frac{\mathrm{\partial }\xi }{\mathrm{\partial }y}& \frac{\mathrm{\partial }\xi }{\mathrm{\partial }z}\\ \frac{\mathrm{\partial }\eta }{\mathrm{\partial }x}& \frac{\mathrm{\partial }\eta }{\mathrm{\partial }y}& \frac{\mathrm{\partial }\eta }{\mathrm{\partial }z}\\ \frac{\mathrm{\partial }\zeta }{\mathrm{\partial }x}& \frac{\mathrm{\partial }\zeta }{\mathrm{\partial }y}& \frac{\mathrm{\partial }}{}\end{array}\end{array}$