Latest revision as of 16:05, 5 February 2020

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

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

Mathematical notation

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

General variables

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

Coordinate systems

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

Kinematics

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

Elasticity

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

Finite element method

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

Modal analysis

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

Numeric integration and SVD

$\beta$ Explicit character parameter
$\mathbf {U}$ Matrix of left side vectors
$\mathbf {V}$ Matrix of right side vectors
${\boldsymbol {\Sigma }}$ Matrix of singular values
$\gamma$ Numerical viscosity
$\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:

(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:

\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:

(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:

(1.4)

Separating real from imaginary parts:

\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 $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 $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:

(1.6)

And the equation of motion results in

(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

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

(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:

(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:

(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

{\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:

(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

(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 ²

(¹)

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

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

$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:

{\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:

(2.2)

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

(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:

{\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}$ :

(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:

(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:

(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:

(2.8)

(2.9)

{\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:

(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:

(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:

(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:

(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:

(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:

(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:

(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:

\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:

(2.20)

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

\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:

(2.22)

An analogous expression can be deduced regarding the elasticity tensor:

(2.23)

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

\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

(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

(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:

(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}$ :

(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:

(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:

{\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:

(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}}}$ .

(3.3)

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

(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}}$ :

\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:

(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

${\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

(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}}}$ :

\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:

(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

{\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

${\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

\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.

+\;{\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} }\;-

\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:

(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:

(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:

(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:

(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:

\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:

{\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:

{\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:

{\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:

(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:

{\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:

{\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} }$ :

\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:

(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:

{\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{matrix}\left[{\begin{array}{l}{dx}\\{dy}\\{dz}\end{array}}\right]=\left[{\begin{array}{lll}{\displaystyle {\frac {\partial x}{\partial \xi }}}&{\displaystyle {\frac {\partial x}{\partial \eta }}}&{\displaystyle {\frac {\partial x}{\partial \zeta }}}\\{\displaystyle {\frac {\partial y}{\partial \xi }}}&{\displaystyle {\frac {\partial y}{\partial \eta }}}&{\displaystyle {\frac {\partial y}{\partial \zeta }}}\\{\displaystyle {\frac {\partial z}{\partial \xi }}}&{\displaystyle {\frac {\partial z}{\partial \eta }}}&{\displaystyle {\frac {\partial z}{\partial \zeta }}}\end{array}}\right]\left[{\begin{array}{l}{d\xi }\\{d\eta }\\{d\zeta }\end{array}}\right]\quad \Rightarrow \quad d{\boldsymbol {x}}=\mathbf {J} ^{e}d{\boldsymbol {\xi }}\;\quad \\\\\left[{\begin{array}{l}{d\xi }\\{d\eta }\\{d\zeta }\end{array}}\right]=\left[{\begin{array}{lll}{\displaystyle {\frac {\partial \xi }{\partial x}}}&{\displaystyle {\frac {\partial \xi }{\partial y}}}&{\displaystyle {\frac {\partial \xi }{\partial z}}}\\{\displaystyle {\frac {\partial \eta }{\partial x}}}&{\displaystyle {\frac {\partial \eta }{\partial y}}}&{\displaystyle {\frac {\partial \eta }{\partial z}}}\\{\displaystyle {\frac {\partial \zeta }{\partial x}}}&{\displaystyle {\frac {\partial \zeta }{\partial y}}}&{\displaystyle {\frac {\partial \zeta }{\partial z}}}\end{array}}\right]\left[{\begin{array}{l}{dx}\\{dy}\\{dz}\end{array}}\right]\quad \Rightarrow \quad d{\boldsymbol {\xi }}={\mathbf {J} ^{e}}^{-1}d{\boldsymbol {x}}\end{matrix}}

(3.28)

Taking a closer look to the elemental Jacobian matrix, one founds that it can be expressed as

\mathbf {J} ^{e}=\left[{\begin{array}{l}{x}\\{y}\\{z}\end{array}}\right]{\begin{bmatrix}{\frac {\partial }{\partial \xi }}&\displaystyle {\frac {\partial }{\partial \eta }}&\displaystyle {\frac {\partial }{\partial \zeta }}\end{bmatrix}}={\boldsymbol {x}}\;{\tilde {\boldsymbol {\nabla }}}_{\boldsymbol {\xi }}^{T}=(\mathbf {X} ^{e}\;{{\mathcal {}}{N}_{(\xi )}^{e}}^{T}){\tilde {\boldsymbol {\nabla }}}_{\boldsymbol {\xi }}^{T}=\mathbf {X} ^{e}\;({\tilde {\boldsymbol {\nabla }}}_{\boldsymbol {\xi }}^{T}{{\mathcal {}}{N}_{(\xi )}^{e}}^{T})=\mathbf {X} ^{e}\;{\tilde {\mathbf {B} }}_{\xi }^{e^{T}}

(3.29)

Where ${\textstyle {\tilde {\boldsymbol {\nabla }}}_{\boldsymbol {\xi }}}$ is the gradient operator for scalar-valued functions in element coordinates and ${\textstyle {\tilde {\mathbf {B} }}_{\xi }^{e}}$ is the matrix of the gradient of ${\textstyle {\mathcal {}}{N}^{e}}$ for scalar-valued functions in the parent domain. The gradient operator for scalar-valued functions in both parent and physical domains ( ${\textstyle {\tilde {\boldsymbol {\nabla }}}_{\boldsymbol {\xi }}}$ and ${\textstyle {\tilde {\boldsymbol {\nabla }}}}$ ) are linked through the transpose of the inverse of the Jacobian matrix:

(3.30)

In the above expression, ${\textstyle {\tilde {\mathbf {B} }}^{e}\in \mathbb {R} ^{n_{d}\times n^{e}}}$ is the matrix of gradient of shape functions for scalar valued functions in global coordinates. In a 3D problem:

{\tilde {\mathbf {B} }}^{e}={\begin{bmatrix}\displaystyle {\frac {\partial N_{1}^{e}}{\partial x}}&\displaystyle {\frac {\partial N_{2}^{e}}{\partial x}}&\cdots &\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial x}}\\\\\displaystyle {\frac {\partial N_{1}^{e}}{\partial y}}&\displaystyle {\frac {\partial N_{2}^{e}}{\partial y}}&\cdots &\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial y}}\\\\\displaystyle {\frac {\partial N_{1}^{e}}{\partial z}}&\displaystyle {\frac {\partial N_{2}^{e}}{\partial z}}&\cdots &\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial z}}\end{bmatrix}}=\left[{\begin{array}{cccc}{{\tilde {\mathbf {B} }}_{1}^{e}}&{{\tilde {\mathbf {B} }}_{2}^{e}}&{\boldsymbol {\cdots }}&{{\tilde {\mathbf {B} }}_{n^{e}}^{e}}\end{array}}\right]

(3.31)

However, the displacement field to be computed is a vector field, so a matrix of symmetric gradient of shape functions for vector-valued functions in the physical domain ${\textstyle \mathbf {B} ^{e}}$ is to be defined. For a 3D problem:

{\mathbf {B} }^{e}={\begin{bmatrix}\displaystyle {\frac {\partial N_{1}^{e}}{\partial x}}&0&0&\cdots &\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial x}}&0&0\\\\0&\displaystyle {\frac {\partial N_{1}^{e}}{\partial y}}&0&\cdots &0&\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial y}}&0\\\\0&0&\displaystyle {\frac {\partial N_{1}^{e}}{\partial z}}&\cdots &0&0&\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial z}}\\\\0&\displaystyle {\frac {\partial N_{1}^{e}}{\partial z}}&\displaystyle {\frac {\partial N_{1}^{e}}{\partial y}}&\cdots &0&\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial z}}&\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial y}}\\\\\displaystyle {\frac {\partial N_{1}^{e}}{\partial z}}&0&\displaystyle {\frac {\partial N_{1}^{e}}{\partial x}}&\cdots &\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial z}}&0&\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial x}}\\\\\displaystyle {\frac {\partial N_{1}^{e}}{\partial y}}&\displaystyle {\frac {\partial N_{1}^{e}}{\partial x}}&0&\cdots &\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial y}}&\displaystyle {\frac {\partial N_{n_{e}}^{e}}{\partial x}}&0\end{bmatrix}}=\left[{\begin{array}{cccc}{{\mathbf {B} }_{1}^{e}}&{{\mathbf {B} }_{2}^{e}}&{\boldsymbol {\cdots }}&{{\mathbf {B} }_{n^{e}}^{e}}\end{array}}\right]

(3.32)

In order to illustrate this section, which may seem a little bit abstract to the reader, let's return to the example presented in Figure 3.1. Imagine a 1D problem, where the domain has been split into smaller straight lines (elements). Each element has at least two nodes where the displacements have to be computed. It is not the only possible case: one can compute the displacements at more than two points by increasing the number of nodes of the element. In 1D, an element with two nodes is called linear; if it has three nodes, quadratic; if four, cubic; and so on. These names refer to the order of the polynomial that has to be used as shape function in each case.

The key is that no matter the length of the element ${\textstyle h^{e}}$ in the physical domain. When mapped to the parent domain, the element boundaries go from ${\textstyle \xi =-1}$ to ${\textstyle \xi =+1}$ . Is in this parent domain where the shape functions ${\textstyle \mathbf {N} ^{e}}$ are defined: as the element length is normalised, the shape functions can be the same for all elements. Recall that in the FEM, a shape function ${\textstyle N_{a}^{e}}$ is restricted to have a value of ${\textstyle \mathbf {1} }$ at node ${\textstyle \mathbf {a} }$ , and ${\textstyle \mathbf {0} }$ at every other node. Notice that in order to accomplish the latter, there will be as many shape functions as nodes. With all these properties in mind, the Figures 3.2 (a), (b) y (c) represent shape functions of different order for a 1D element defined in the parent domain.


(a)	(b)

(c)
Figure 3.2: 1D shape functions for linear (a), quadratic (b) and cubic (c) isoparametric elements.

The same reasoning that has been followed for a 1D element can be developed for 2D and 3D elements. In those cases, many different types of elements can be defined, depending on their geometric properties and number of nodes. The shape functions will then have two or three spatial variables, and are to be defined for each type of element. For instance, in order to complete the previous example, let's consider the case of a linear triangular element (Figure 3.3).

Figure 3.3: Representation of a triangular linear element in both parent and physical domains

In the parent domain, the nodes are always listed in counter-clockwise direction. As can be seen, no matter the shape of the triangle in the physical domain, once it is mapped into the parent space all the elements are geometrically identical. The matrix of shape functions of the above element and its gradient in the parent domain are, then:

(3.33)

3.3.1 The Gauss Quadrature integration method

Although the domain has been divided into smaller elements, the integrals presented earlier are still difficult to compute, as those elements have a wide variety of different boundaries. If the geometry of the element is sort of complex (nothing unusual in 3D), the integrals would have to be evaluated along those boundaries, which also change from one element to another. The solution to this problem has been presented in the previous section: if the element is mapped into the parent domain, its geometric properties are kind of normalised. The integrals can then be evaluated in the parent domain, and later on be transformed into global coordinates by using the Jacobian matrix (recall equation 3.30).

Now the question lies in which method use to evaluate those integrals. For sure, the analytic solution is discarded as it would take too much time. From among all the numeric integration methods, one of the most widespread is Gauss Quadrature. It is a fast and simple method, that give great approximations with few operations. Basically, it reduces the integral to a summation of function values at specific points (Gauss points), multiplied by a certain factor called weight. The great advantage of this method is that as all the elements of the same type are mapped equally in the parent domain, the values of the weights and gauss points will be the same for all elements.

Imagine that the integral to evaluate ${\textstyle I}$ is of some function ${\textstyle f:{\overline {\mho }}^{e}\rightarrow \mathbb {R} }$ over the physical domain. First, the integral has to be mapped into the parent domain through a change of variable, accomplished by using the Jacobian matrix of the transformation ${\textstyle {\boldsymbol {x}}={\boldsymbol {x}}(\xi )}$ :

(3.34)

Where ${\textstyle J}$ stands for the determinant of the Jacobian. In 3D, ${\textstyle \mathrm {d} \Omega _{\xi }=\mathrm {d} \xi \;\mathrm {d} \eta \;\mathrm {d} \zeta }$ . As Gaussian rules for integrals in several dimension are constructed by using 1D Gaussian rules on each coordinate separately, hereinafter only the one-dimensional Gauss Quadrature will be considered. The goal of the method is to find a set of Gauss points ${\textstyle \xi _{g}}$ and weights ${\textstyle w_{g}}$ such that:

(3.35)

The above expression give exact results only if ${\textstyle J^{e}(\xi )f(\xi )}$ is a polynomial. If not, only approximate results are obtained, whose accuracy improves as more Gauss points are used. Supposing that the integrand ${\textstyle q(\xi )=\mathbf {J} ^{e}(\xi )f(\xi )}$ is a polynomial of order ${\textstyle p}$ :

(3.36)

The values of the Gauss points and weights are obtained by exactly evaluating each monomial ${\textstyle \xi ^{i}}$ , with ${\textstyle i=1,2\ldots p}$ :

\sum _{g=1}^{m}w_{g}\xi _{g}^{i}=\left.{\frac {\xi ^{i+1}}{i+1}}\right|_{\Omega _{\xi }}^{}\;{\xrightarrow {=}}\;\;=\;\left\{{\begin{array}{ll}{0}&{{if}\quad i=1,3,5\ldots }\\{\frac {2}{i+1}}&{{if}\quad i=0,2,4,6\ldots }\end{array}}\right.

(3.37)

The previous expressions lead to a linear system of ${\textstyle p+1}$ equations with ${\textstyle 2m}$ unknowns. If ${\textstyle q}$ results to be a polynomial, the exact integration requires from ${\textstyle m\geq {\frac {p+1}{2}}}$ points. Gauss points and weights will be different depending on the number of points ${\textstyle m}$ used. The values of the Gauss points and weights up to 5 points are presented in Table 3.1.

Table 3.1. Gauss points positions and weights [9]

4. FEM programming

Now that the elastodynamic model for rotating structures has been formulated in terms of the FEM, it is time to tackle its programming. This section covers the different processes involved in the simulation and how the ${\textstyle MATLAB}$ solver works. In this chapter, only the static case is to be addressed, as the methods used to solve the dynamic problem will be reviewed later in the report.

4.1 Preprocessing

All the finite element simulations, carried out or not by commercial codes, start with the model definition, a critical phase widely known as preprocessing. This first step begins with the definition of the geometric domain of the problem: in this case, the structure. Nowadays, the most widespread method to define the geometry is by using computer-aided design (CAD) software. Examples of commercial CAD programs are CATIA, Rhino and SolidWorks.

Once the geometry is defined, it is exported to what is known as mesh generation software. The final goal of the latter is to create a mesh: a subdivision of the geometry into discrete cells known as elements. As it has been commented in the previous chapter, many different types of elements exist, so its type, length, area and other properties are to be defined. Notice that meshing is not a trivial process as geometries can be quite intricate. In order to accurately capture the input domain geometry, the mesh generation is often aided by computer algorithms. In this step, the domain boundary is also to be defined. Once the mesh has been generated, the software computes what is known as connectivity matrix, which contains all the information about the mesh and defines the relations between elements and nodes (see for instance the numbers printed in Figure 4.1a).

The connectivity matrix ${\textstyle \mathbf {CN} \in \mathbb {R} ^{n_{el}\times n^{e}}}$ , also known as element connectivity array, stores all the element definitions. It is an array of integers numbers where the ${\textstyle \mathbf {(e\,,\,a)} }$ cell contains the global number of the ${\textstyle \mathbf {a} |th}$ node of element ${\textstyle \mathbf {e} }$ . It is an essential matrix as it stores all the information regarding the mesh, and it is crucial when transforming from global to local coordinates.

Quite often, commercial programs allow both designing and meshing. This is in fact the case of GiD, the pre and postprocessing tool used during this project. However, its designing tools are quite basic, and so it admits the importation of models from other CAD softwares. Nonetheless, those tools are quite useful when editing an imported geometry in order to simplify the meshing procedure. Figure 4.1 illustrates the meshing of a 2D figure with GiD, using triangular and quadrilateral elements. In some cases, the mesh can be very simple (for instance a grid in Figure 4.1a), whereas other geometries require from more complex meshes.

The next step is to define, for each element, material properties (density, elastic coefficients, heat conductivity, etc.) and loadings (volumetric forces, internal heat sources, etc.). The final part of the preprocessing phase consists in defining boundary conditions, which in our case is to prescribe displacements and tractions along the boundary. Commercial softwares usually include a wide library of material properties and loading scenarios to carry out these last steps.


(a)	(b)
Figure 4.1: Mesh representation using quadrilateral (a) and triangular (b) elements. In the case of the circle, the mesh has been generated using in-built algorithms.

As GiD does not implement this capability, this last part of the preprocessing is carried out with ${\textstyle MATLAB}$ . Once the mesh is created, GiD exports a raw file (.msh extension) including the connectivity matrix, the nodal coordinates, the element typology and boundary information. Those are later recovered and imported to ${\textstyle MATLAB}$ . Material properties are manually assigned to each element, as well as the nodal values of body forces and prescribed boundary displacements and tractions.

4.2 Solver

Once the formulation of the problem is completed and all the information regarding structure's geometry, material properties, loads and boundary conditions is gathered in a ${\textstyle MATLAB}$ file, it is time to solve the elastostatic equations. In fact, the equations to be solved correspond to the approximated weak form of the problem stated in terms of the FEM (see equation 3.21). In this chapter, only the static form of the problem is to be tackled — that is, ${\textstyle {\ddot {\mathbf {d} }}={\dot {\mathbf {d} }}=0}$ :

(4.1)

It is a first step in the understanding of the whole phenomena, and the developed code is going to be re-used later on when studying the dynamic case. The code is developed in ${\textstyle MATLAB}$ , with its main working algorithm schematised in Figure 4.2.

Figure 4.2:

main code scheme for elastostatic problems

As it can be suspected from the previous figure, the main task of the solver is to compute the matrix ${\textstyle \mathbf {K} }$ and vector ${\textstyle \mathbf {F} }$ from equation 3.13. Using the local formulation presented in the previous chapter, the computation of the latter can be split into a simpler problem, regarding only elemental matrices ${\textstyle \mathbf {K} ^{e}}$ and vectors ${\textstyle \mathbf {F} ^{e}}$ . Once they are obtained for all the elements of the domain, the elemental matrices are joined together into a global stiffness matrix ${\textstyle \mathbf {K} }$ and force vector ${\textstyle \mathbf {F} }$ through a process called assembly. When ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {F} }$ are known, equation 4.1 results in a linear system of algebraic equations, easily handled with native ${\textstyle MATLAB}$ functions. A further explanation regarding elemental matrix computation and assembly can be found in Annex B.

Once the displacement field is known, the final step of the solver consists in finding the reactions and the strain field. Computation of elastic reactions is straightforward, as they are obtained from equation 3.21. On the other hand, the strain field is to be computed in the local frame, element by element. Strains are defined as the spatial derivative of displacements (see equation A.7), which can be written in terms of the FEM as:

(4.2)

Where ${\textstyle \mathbf {A} }$ represents the assembly operator. As ${\textstyle \mathbf {B} ^{e}}$ is approximated through Gauss Quadrature, strains and stresses are only computed at the Gauss points ${\textstyle \xi _{g}}$ .

4.2.1 Efficient assembly

In the FEM, computation of global matrices is never done in a single step. Instead, elemental matrices are computed element by element, and then joined together through a process called assembly, constructing this way the global matrix. The intuitive method to perform the assembly is quite simple to program: elemental matrices are computed one by one, and then allocated into a bigger, global matrix using the connectivity array. The connectivity array ${\textstyle \mathbf {CN} }$ links the local numeration of the degrees of freedom with the global one, specifying the row and column of the global matrix where each component of the elemental matrix has to be placed. For a one-dimensional case, the classical assembly of a matrix ${\textstyle \mathbf {M} }$ can be programmed as follows:

for e = 1 : n_el %Loop over all the elements
   M_e = Compute_Me(. . .); %Function that computes the elemental matrix
   for a = 1 : n_e %Loop over the nodes of the element e
       for b = 1 : n_e
          A = CN(e,a); B = CN(e,b); %Transformation from local to global
          M_glo(A,B) = M_glo(A,B)+M_e(a,b); %Allocate M_e into the global matrix
       end
   end
end

The previous algorithm has two main drawbacks. The first one is related with the memory large matrices take up. Each cell of a ${\textstyle MATLAB}$ matrix is defined as double, a type of variable that uses 8 bytes of memory. Using the traditional approach, a 3D system with 30000 nodes would involve matrices with a size bigger than 7.5 gigabytes, which is clearly infeasible. But if one takes closer look to the matrices involved in the FEM, will discover that most of the cells are filled with zeros. Matrices in which most of the elements are zero are known as sparse matrices Figure 4.3.

Figure 4.3: Representation of a sparse matrix. Blue squares are non-zero values.

For sure, storing the zeros of a sparse matrix as double variables is a waste of memory. That is why most of the programming languages have in-built algorithms to store sparse matrices in more efficient ways. The key is to store only the non-zero values and their position in the matrix. Using this approach, systems with a large number of degrees of freedom can be tackled without memory problems.

The second disadvantage of the classical assembly method is related with its computational cost. A closer look to the presented algorithm reveals three nested loops. A local matrix has to be computed for each element, and then placed into the global matrix cell by cell. In ${\textstyle MATLAB}$ , this iteration is rather slow. For the static case, where matrices are computed once, it is not a big deal. But when dealing with the dynamic case, global matrices are to be computed repeatedly, increasing considerably the running time of the solver. This problem is solved by using a vectorised version of the assembly algorithm.

The key of the vectorised algorithm is that computes and assembles all the elemental matrices at once, avoiding iterations [10]. It is more difficult to program than the classical algorithm, but its advantages are considerable. For instance, Figure 4.4 reveals an outstanding decrease of the computation time for 2D systems. As it can be spotted, although for few degrees of freedom the classical algorithm is quicker, the vectorised algorithm proves to be the most efficient solution as the size of the system grows.

The vectorised assembly algorithm has been adapted for the rotating case from a set of codes given by prof. Hernández J., which were used in the context of composite materials. An exemplification of the vectorised assembly for a simple mesh can be found in Annex C.

Figure 4.4: Comparison of the computation time of

,

and

using the classical and vectorised assembly (logarithmic scale).

4.2.2 Approximated one-dimensional stationary solution

To conclude the solver section, let's exemplify the explained concepts approximating the solution of the system presented in 2.26. It is the simplest case to tackle: a 1D blade rotating at constant speed with no vibration at all. It can be interpreted as a static elastic problem in the rotational frame, where deformations can only occur in the axial direction ${\textstyle \varsigma }$ . The goal is to find the approximated solution using the FEM, and assess its validity by comparing it with the analytic result obtained in 2.5.2. Recall the differential equation to solve:

(4.3)

Using the FEM formulation, an approximated weak form of the equation can be written:

(4.4)

By replacing the test and trial function with its approximations in terms of the shape functions, the above equation yields to:

{\overset {\mathbf {K} }{\overbrace {\left({\underset {\mathbf {K} _{\hbox{static}}}{\underbrace {\int _{\mho }\mathbf {B} ^{T}\,\,{\textrm {E}}\;\;\mathbf {B} \;\;d\varsigma } }}-{\underset {\mathbf {K} _{\hbox{rot}}}{\underbrace {\int _{\mho }\mathbf {N} ^{T}\;\rho \;\Omega ^{2}\;\mathbf {N} \;\;d\varsigma } }}\right)} }}\;\mathbf {d} ={\overset {\mathbf {F} \;=\;\mathbf {F} _{\hbox{rot}}}{\overbrace {\int _{\mho }{\mathbf {N} }^{T}\;\rho \;\Omega ^{2}\;\;\varsigma \;d\varsigma } }}

(4.5)

Taking into account that displacements are prescribed along the set ${\textstyle \mathbf {r} }$ and unknown in ${\textstyle \mathbf {l} }$ , the above equation can be split into block matrices:

(4.6)

Previous to the solving of the above system of algebraic linear equations, global matrices ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {F} }$ are to be known. As it has been discussed, the best choice is to decompose the latter into elemental matrices:

(4.7)

Elemental matrices are to be calculated using the approach reviewed in Annex B. The first step is to define the type of element in the parent domain. For simplicity, let's assume linear elements. Recall the shape functions in the parent domain for a linear 1D element from Figure 3.2(a):

(4.8)

To compute ${\textstyle \mathbf {B} _{\xi }^{e}}$ , the Jacobian of the transformation ${\textstyle {\frac {\mathrm {d} \xi }{\mathrm {d} \varsigma }}}$ is to be found:

(4.9)

Where ${\textstyle h^{e}}$ represents the length of the element. Then:

(4.10)

With all, the computation of the elemental stiffness matrix is reduced to:

\mathbf {K} _{\hbox{static}}^{e}=\int _{\mho }\mathbf {B} ^{e^{T}}\,\,{\textrm {E}}\;\;\mathbf {B} ^{e}\;\;d\varsigma =\int _{-1}^{1}{\frac {1}{h^{e}}}\left[-1\;\;\;\;+1\right]^{T}{\textrm {E}}\;\;{\frac {1}{h^{e}}}\left[-1\;\;\;\;+1\right]\;{\frac {h^{e}}{2}}d\xi ={\frac {\textrm {E}}{h^{e}}}\;{\begin{bmatrix}1&-1\\-1&1\end{bmatrix}}

(4.11)

Recall that for the calculation of the elemental stiffness matrix, integrals are evaluated exactly. This is not the case of the force vector, which is approximated trough Gauss Quadrature. Replacing ${\textstyle \rho \;\Omega ^{2}\;\varsigma }$ by ${\textstyle \mathbf {f} }$ :

={\frac {h^{e}}{8}}\sum _{g=1}^{m}w_{g}\rho ^{e}\;\Omega ^{2}\;{\begin{bmatrix}(1-\xi _{g})^{2}&(1-\xi _{g})\;(1+\xi _{g})\\(1-\xi _{g})\;(1+\xi _{g})&(1+\xi _{g})^{2}\end{bmatrix}}{\begin{bmatrix}\varsigma _{1}^{e}\\\varsigma _{2}^{e}\end{bmatrix}}

(4.12)

As ${\textstyle \mathbf {N} ^{e}}$ is linear in ${\textstyle \xi }$ and ${\textstyle \mathbf {f} }$ is also linear in ${\textstyle \varsigma }$ , three Gauss points ( ${\textstyle m=3}$ ) will suffice to evaluate the above integral exactly. The values of Gauss points and weights can be found in table 3.1. Once the elemental stiffness matrix and force vector are known, they are computed for each element of the domain and assembled into global matrix ${\textstyle \mathbf {K} }$ and vector ${\textstyle \mathbf {F} }$ . The last step before finding displacement field is to solve the system of equations stated in 4.6.

In order to compare the results with the exact solution, let's consider the same conditions as in the analytic case. The structure of study will be a blade with constant cross-section and ${\textstyle 4\;m}$ radius, rotating at ${\textstyle 60\;{\hbox{rad/s}}}$ . For an aluminium bar ( ${\textstyle E=70\;GPa}$ , ${\textstyle \rho =2700\;kg/m^{3}}$ ), the results are displayed in Figure 4.5 and compared with the exact solution. Recall that computation of stresses once nodal displacements are known is achieved using expression 4.2.

Figure 4.5: Approximated solution of equation 2.26 using 70 elements, compared with the exact result.

In the FEM, another input parameter to define is the number of elements into which the domain is divided. Theoretically, the approximated solution converges into the exact one as the number of elements increases. In order to assess this convergence, let's define the error associated to the approximated solution ${\textstyle \mathbf {e} }$ using the L2 norm of the distance between ${\textstyle \mathbf {u} _{\varsigma }}$ and ${\textstyle \mathbf {u} _{\varsigma }^{h}}$ :

(4.13)

In a similar way, one can define the error ${\textstyle \mathbf {e'} }$ associated to the derivative of ${\textstyle \mathbf {u} }$ , which in this case is linked with stresses. In order to evaluate those integrals, the domain is again split into elements where Gauss Quadrature is applied. As quite often the exact solution is not a polynomial, high order Gauss rules are needed. In general, the logarithm of the error varies linearly with the element size, and its slope depends on the order of the element:

(4.14)

Where ${\textstyle \alpha }$ is known as the rate of convergence of the element, while ${\textstyle C}$ is an arbitrary constant and ${\textstyle h}$ the element length (assumed constant). The parameter ${\textstyle \alpha }$ is equal to ${\textstyle p+1}$ , where ${\textstyle p}$ is the order of the used shape function. For the derivative of the error, the slope of convergence is one order lower. Recall that if obtaining the analytic exact solution were impossible, the computation of the error in the way that has been developed would not be practicable. In that case, in order to check the convergence of the method, relative error should be computed. That is, to compare the error of the function obtained using ${\textstyle n_{el}}$ elements, with a previous function obtained with ${\textstyle n_{el}-k}$ elements.

4.3 Postprocessing

The last step of every FEM simulation is analysis and evaluation of the solution results, often referred as postprocessing. Once the solution of the problem is computed, it is exported to a postprocess program that interprets, prints and plots the results. Among many other operations, these programs are used to draw the deformed structural shape, produce colour-coded plots (representing displacements, strains and stresses) and animate the dynamic behaviour of the system, if that were the case. Equilibrium check (assessing if reactions counteract external forces, for example) is also part of the postprocess stage.

In the previous example, results were graphed in a simple plot. This approach is unattainable for multidimensional problems with non-trivial geometry, and thus external programs are to be used. In this project, the postprocessing is carried out by GiD, the same program used for the mesh definition. To illustrate GiD capabilities, let's consider a fixed-free cantilever beam under the effect of a distributed load ${\textstyle w=-20\;KN/m}$ along its upper side. The beam is shaped as a long rectangle ( ${\textstyle 0,2\times 2m}$ ) and meshed using 1000 linear quadrilateral elements. Figure 4.6 represents its deformed configuration, and stresses along the beam are represented using colours:

Figure 4.6: Deformed shape (

) and stress (

) colour map of a 2D beam under the effect of a distributed load on its upper side.

The colour legend has not been displayed as it is not from interest in this example. Remember that although displacement is continuous along the domain, its derivatives (stress and strains) are only computed at Gauss points. Modern postprocessing programs have in-built tools that blur the latter in order to hide discontinuities. The engineer has to be aware that misuse of blurring tools can induce notable errors.

To both asses the accuracy and validate the 2D solver, FEM results can be contrasted with the analytic solution obtained from beam theory. In particular, vertical displacement along the mean line of the beam is the variable that is being compared in Figure 4.7. Analytic results are those of a beam with squared cross section, whose deflection can be expressed as:

(4.15)

Figure 4.7: Comparison of the approximated and theoretical vertical displacement of a beam under the effect of a distributed load.

5. Dynamic response

So far, the problem regarding static deformations in the noninertial frame is fully settled. However, although it can be used to determine initial conditions, the developed model is unable to account for time-dependant variables such as velocity or acceleration in the rotating frame. In other words: it cannot capture the phenomena of vibration. In order to study the dynamic response of the system, new approaches are to be considered. Hence, this chapter will review two main methodologies to tackle time-dependant problems: modal decomposition analysis and numerical time integration. A brief description regarding the singular value decomposition (SVD) will be presented at the end of this chapter, presenting a promising method to handle nonlinear dynamic systems.

5.1 Modal decomposition analysis

Modal analysis is one of the most widespread techniques to tackle structural vibration problems, for both discrete and continuous systems. The key of this method is that it studies the problem in the frequency domain, rather than in the temporal one. To start with, recall the equation describing a rotating structure in the noninertial frame of reference:

\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} }})

(5.1)

Where ${\textstyle \mathbf {K} }$ , ${\textstyle \mathbf {D} }$ and ${\textstyle \mathbf {F} }$ are the sum of two contributions: a first one typical of static conditions, and a particular second one accounting for the effect of rotation. For the sake of notational simplicity, on now on let's assume the following redefinitions:

{\begin{matrix}\mathbf {M} _{\mathbf {ll} }\rightarrow \mathbf {M} \;\;,\;\;\;\;\;\;\mathbf {K} _{\mathbf {ll} }\rightarrow \mathbf {K} \;\;,\;\;\;\;\;\;\mathbf {D} _{\mathbf {ll} }\;\rightarrow \mathbf {D} \;\;,\;\;\;\;\;\;\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} }})\rightarrow \mathbf {F} \\\\\mathbf {d} _{\mathbf {l} }\rightarrow \mathbf {d} \;\;,\;\;\;\;\;\;\mathbf {d} _{\mathbf {l} }^{0}\rightarrow \mathbf {d} ^{0}\;\;,\;\;\;\;\;\;{\dot {\mathbf {d} }}_{\mathbf {l} }^{0}\rightarrow {\dot {\mathbf {d} }}^{0}\end{matrix}}

(5.2)

The advantage of the presented reformulation is that allows equation 5.1 to be expressed as:

(5.3)

With ${\textstyle \mathbf {d} :[0,{\textrm {T}}]\rightarrow \mathbb {R} ^{n}}$ and subject to initial conditions ${\textstyle \mathbf {d} _{(0)}=\mathbf {d} ^{0}}$ and ${\textstyle {\dot {\mathbf {d} }}_{(0)}={\dot {\mathbf {d} }}^{0}}$ . The previous expression can be thought as a second order ordinary linear differential equation. From calculus, it is known that the solution of this kind of ODE is the sum of a particular ${\textstyle \mathbf {d} ^{\mathbf {p} }}$ and homogeneous ${\textstyle \mathbf {d} ^{\mathbf {z} }}$ solution:

(5.4)

The homogeneous solution of the system is found by making ${\textstyle \mathbf {F} =0}$ , whereas the particular solution is to be computed for the specific case of study. One of the main drawbacks of the modal analysis is that, as it works in the frequency domain, forcing functions need to be harmonic: ${\textstyle \mathbf {F} ={\overline {\textrm {F}}}\,\sin(\,{\textrm {wt}}\,)}$ . This is not the case of the system of study, where forces due to rotation are either constant or increasing in time, and external forces can adopt any arbitrary shape. One way of tackling this problem is to decompose ${\textstyle \mathbf {F} }$ into the sum of harmonic functions through the Fourier transform. As the system of study is linear, its response can be expressed as the sum of the responses due to each of the obtained harmonic forcing functions.

If ${\textstyle \mathbf {F} }$ is difficult to pose in terms of harmonic functions, another approach consists in splitting ${\textstyle \mathbf {F} }$ into the sum of ${\textstyle n}$ impulses. The response of the system to an impulse is identical to the free response with certain initial conditions [11]. Again, the system is studied under the effect of each of those individual impulses, and taking advantage of the linearity of the system, the solution is expressed as the sum of the latter individual responses. As modal study of structures under a generic forcing function ${\textstyle \mathbf {F} }$ can easily become complicated, the present section will only cover the homogeneous solution of equation 5.3.

5.1.1 Undamped free vibrations

In order to introduce the basics of the modal analysis, let's consider the simplest case to handle: a system where both ${\textstyle \mathbf {F} }$ and ${\textstyle \mathbf {D} }$ are zero, so damping and external forces play no role whatsoever. With this assumption, equation 5.3 reduces to:

(5.5)

The modal analysis starts supposing a solution with the following form:

(5.6)

Where ${\textstyle {\textrm {'''q'''}}_{(t)}}$ is a known solution for homogeneous second order differential equations. Since ${\textstyle {\ddot {\mathbf {q} }}=-{\boldsymbol {\omega }}^{2}\;{\mathbf {q} }}$ , the above equation can be written as:

(5.7)

The preceding equation is known as the eigenvalue problem. The trivial solution ${\textstyle {\boldsymbol {\Phi }}=0}$ is not from interest, whereas the non-trivial solution is given by a set of eigenvalues ${\textstyle {\boldsymbol {\lambda }}={\boldsymbol {\omega }}^{2}}$ , where ${\textstyle {\boldsymbol {\omega }}}$ represents the diagonal matrix of natural or resonant vibration frequencies. For a system with ${\textstyle \mathbf {n} }$ degrees of freedom, ${\textstyle \mathbf {n} }$ eigenvalues exist that satisfy the above expression. Being ${\textstyle j}$ an index that goes from ${\textstyle 1}$ to ${\textstyle \mathbf {n} }$ , the eigenvalue problem can be interpreted as a system of equations of the form

(5.8)

Where ${\textstyle {\boldsymbol {\Phi }}_{j}}$ is referred as the eigenvector associated to the eigenvalue ${\textstyle {\lambda }_{j}}$ . Notice that eigenvectors do never define absolute values of the different independent coordinates of the system, but inform about the phase and amplitude relations between degrees of freedom. In structural problems, eigenvectors are called natural modes or mode shapes, and can be though as possible deformed configurations oscillating at a given frequency stated by its corresponding eigenvalue. In matrix form, eigenvalues and eigenvectors are expressed as:

{\boldsymbol {\lambda }}={\begin{pmatrix}\omega _{1}\,^{2}&0&\cdots &0\\0&\omega _{2}\,^{2}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &\omega _{n}\,^{2}\end{pmatrix}}\;\;\;,\;\;\;\;\;\;{\boldsymbol {\Phi }}={\begin{bmatrix}{\boldsymbol {\Phi }}_{1}&{\boldsymbol {\Phi }}_{2}&\cdots &{\boldsymbol {\Phi }}_{n}\end{bmatrix}}={\begin{pmatrix}\Phi _{11}&\Phi _{12}&\cdots &\Phi _{1n}\\\Phi _{21}&\Phi _{22}&\cdots &\Phi _{2n}\\\vdots &\vdots &\ddots &\vdots \\\Phi _{n1}&\Phi _{n2}&\cdots &\Phi _{nn}\end{pmatrix}}

(5.9)

The main advantage of this method is that natural frequencies and mode shapes are relatively easy to obtain. ${\textstyle MATLAB}$ , for instance, has in-built functions that compute a set of ${\textstyle \mathbf {m} \leq \mathbf {n} }$ eigenvectors with the smallest associated frequency. The key of this approach is that there is no need to compute all the ${\textstyle \mathbf {n} }$ eigenvalues and eigenvectors, and instead the solution of the system can be approximated using a small number of them. In structural analysis, shape modes associated to low frequencies are the most dangerous, and the ones from interest when studying vibration. Recall that ${\textstyle {\boldsymbol {\Phi }}}$ and ${\textstyle {\boldsymbol {\lambda }}}$ only depend on the matrices ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {M} }$ , which can be computed using the static FEM program previously presented. For instance, mode shapes associated to the four lowest natural frequencies of the cantilever beam defined in 4.6 are displayed in Figure 5.1.

Figure 5.1: Shape modes associated to the first four lower natural frequencies of a 2D fixed-free beam

An important property of eigenvectors is that they are orthogonal with respect to ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {M} }$ :

(5.10)

The previous equation can only be applied if the matrices involved in the system are symmetric and positive semidefinite, and implies that the matrix of natural modes ${\textstyle {\boldsymbol {\Phi }}}$ diagonalizes both ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {M} }$ . It is a common approach to normalise ${\textstyle {\boldsymbol {\Phi }}}$ such that ${\textstyle {\boldsymbol {\Phi }}_{i}^{T}\mathbf {M} {\boldsymbol {\Phi }}_{i}=\mathbf {I} }$ , which results into ${\textstyle {\boldsymbol {\Phi }}_{i}^{T}\mathbf {K} \,{\boldsymbol {\Phi }}_{i}={\boldsymbol {\omega }}^{2}}$ . Introducing the generalised coordinates ${\textstyle \mathbf {q} \in \mathbb {R} ^{n}}$ such that ${\textstyle \mathbf {d} ={\boldsymbol {\Phi }}\;\mathbf {q} }$ , and pre-multiplying by ${\textstyle {\boldsymbol {\Phi }}_{i}^{T}}$ equation 5.5:

(5.11)

The above expression is, in fact, a system of uncoupled differential equations subject to the following initial conditions:

(5.12)

Recalling the definition of ${\textstyle \mathbf {q} }$ stated in equation 5.6, the exact solution of the system can be expressed in uncoupled form as:

(5.13)

In order to postprocess the results, displacements of every degree of freedom are to be computed at different values of time and joined together into a matrix. The time between two consecutive samples is called time-step, and in order to capture the vibration phenomena it has to be less than half the period of the highest frequency of interest. Once displacements are computed, they can be animated through GID.

\mathbf {d} ={\begin{bmatrix}\mathbf {d} _{({\textrm {t}}_{1})}&\mathbf {d} _{({\textrm {t}}_{2})}&\cdots &\mathbf {d} _{\textrm {(T)}}\end{bmatrix}}={\begin{pmatrix}{{\textrm {d}}_{1}}_{({\textrm {t}}_{1})}&{{\textrm {d}}_{1}}_{({\textrm {t}}_{2})}&\cdots &{{\textrm {d}}_{1}}_{({\textrm {T}})}\\{{\textrm {d}}_{2}}_{({\textrm {t}}_{1})}&{{\textrm {d}}_{2}}_{({\textrm {t}}_{2})}&\cdots &{{\textrm {d}}_{2}}_{({\textrm {T}})}\\\vdots &\vdots &\ddots &\vdots \\{{\textrm {d}}_{n}}_{({\textrm {t}}_{1})}&{{\textrm {d}}_{n}}_{({\textrm {t}}_{2})}&\cdots &{{\textrm {d}}_{\mathit {n}}}_{({\textrm {T}})}\end{pmatrix}}

(5.14)

As it has been pointed before, sometimes it is enough with an approximation of the solution, obtained using only ${\textstyle \mathbf {m} <\mathbf {n} }$ terms. A good strategy to asses which shape modes are the most relevant for a given structure is to use the normalised initial amplitude associated to each mode as guide:

(5.15)

For the presented 2D beam, normalised initial amplitudes are plotted in a bar graph in Figure 5.2 for the six shape modes with lower frequency. It can be seen that the most excited mode shape is the one corresponding to pure bending ( ${\textstyle {\boldsymbol {\Phi }}_{1}}$ ), which makes sense since the beam is under the effect of a distributed load on its upper side. Notice that in order to find this out, initial conditions need to be defined. Different modes can dominate the same structure under distinct initial conditions.

Figure 5.2: Normalised initial amplitude versus mode shape number for a 2D fixed-free beam under the effect of a distributed load

5.1.2 Damped free vibrations

Modal analysis can handle damped systems only if some requirements regarding the damping matrix are met: ${\textstyle \mathbf {D} }$ needs to be either diagonal or a linear combination of ${\textstyle \mathbf {M} }$ and ${\textstyle \mathbf {K} }$ . The last assumption is known as Rayleigh damping, and has already been presented in equation 3.15. The equation describing a system under Rayleigh damping is written as:

(5.16)

Introducing the change of variable ${\textstyle \mathbf {d} ={\boldsymbol {\Phi }}\;\mathbf {q} }$ , and pre-multiplying by ${\textstyle {\boldsymbol {\Phi }}_{i}^{T}}$ :

{\overset {\mathbf {I} }{\overbrace {{\boldsymbol {\Phi }}^{T}\mathbf {M} {\boldsymbol {\Phi }}} }}\,{\ddot {\mathbf {q} }}+({\overline {\alpha }}\;{\overset {\mathbf {I} }{\overbrace {{\boldsymbol {\Phi }}^{T}\mathbf {M} {\boldsymbol {\Phi }}} }}+{\overline {\beta }}\;{\overset {{\boldsymbol {\omega }}^{2}}{\overbrace {{\boldsymbol {\Phi }}^{T}\mathbf {K} \;{\boldsymbol {\Phi }}} }})\,{\dot {\mathbf {q} }}+{\overset {{\boldsymbol {\omega }}^{2}}{\overbrace {{\boldsymbol {\Phi }}^{T}\mathbf {K} \;{\boldsymbol {\Phi }}} }}\,\mathbf {q} \;\;\;\Rightarrow \;\;\;{\ddot {\mathbf {q} }}+\left({\overline {\alpha }}+{\overline {\beta }}\,{\boldsymbol {\omega }}^{2}\right){\dot {\mathbf {q} }}+{\boldsymbol {\omega }}^{2}\mathbf {q} =\mathbf {0}

(5.17)

Which is a system of uncoupled differential equations. Defining the damping ratio associated to the ${\textstyle i|th}$ eigenvalue ${\textstyle {\overline {\xi }}_{i}}$ and the damped oscillation frequency ${\textstyle {\overline {\omega _{i}}}}$ :

(5.18)

The system of equations can be decoupled into:

(5.19)

Up to this point, two different situations may occur depending on the value of ${\textstyle {\overline {\xi }}_{i}}$ . If ${\textstyle 0\leq {\overline {\xi }}_{i}<1}$ , then ${\textstyle {\overline {\omega }}_{i}}$ is real valued and the exact expression for the generalised coordinates is:

\mathbf {q} _{i}=e^{-{\overline {\xi }}_{i}\omega _{i}{\textrm {t}}}\left(\mathbf {q} _{i}^{0}\cos {\overline {\omega }}_{i}{\textrm {t}}+{\frac {{\dot {\mathbf {q} }}_{i}^{0}+{\overline {\xi }}_{i}\omega _{i}\,\mathbf {q} _{i}^{0}}{{\overline {\omega }}_{i}}}\sin {\overline {\omega }}_{i}{\textrm {t}}\right)

(5.20)

This response is known as subcritical oscillation, and can be interpreted as an undamped free vibration whose amplitude decreases exponentially in time. On the other hand, if ${\textstyle {\overline {\xi }}_{i}>1}$ , damped oscillation frequencies are no longer real valued and generalised coordinates are expressed as:

(5.21)

This is a non-oscillating response known as overcritical motion. In both cases, physical displacements can be computed in the same way as in the undamped case:

(5.22)

5.1.3 Modal analysis of rotating structures

The presented method of modal decomposition can be applied to the case of study of rotating structures, but subject to some limitations. The approach that has been developed previously is only valid if the matrices involved in the system are constant. Recall from equation 3.13 that for rotation problems, the stiffness and damping matrices are the addition of a static component, which is constant, and a rotating one that depends on angular velocity and acceleration. The latter varies with time if the angular acceleration is different from zero. Thus, the only possibility for the system to be tackled through modal analysis is to rotate under constant angular velocity ${\textstyle \Omega }$ . In that case, the eigenvalue problem is written as:

(5.23)

And if the problem is either 1D or 2D:

(5.24)

The above expression unfolds the softening effect predicted in 3.17, which produces a decrease in the natural frequencies as angular velocity increases. The relationship between natural frequency of the static and rotating case is:

(5.25)

According to the above equation, lower frequencies are the most affected by the softening effect. This behaviour can be represented through a Campbell diagram (Figure 5.3), where natural frequencies of a rotating structure are plotted as a function of the angular speed. This softening effect, however, is rarely seen in real scenarios: it is a consequence of the limitations of the model, which is not able to account for a counteracting stiffening effect. The decrease of the natural frequencies emerges in fact due to the non-inertial character of the frame of study.

Supposing that the 2D beam studied in Figure 5.1 undergoes a rotation at constant angular velocity, then the natural frequencies linked with the presented mode shapes would decrease as the spinning speed grows in the way it is displayed in the following graph.

Figure 5.3: Campbell diagram representing the natural frequencies of a 2D beam rotating at constant speed

5.1.4 Limitations

Until now, only stationary scenarios have been considered as modal analysis is not capable of accounting the effect of angular acceleration. Another drawback of this method is related with the requirement that ${\textstyle \mathbf {D} }$ needs to be a linear combination of ${\textstyle \mathbf {K} }$ and ${\textstyle \mathbf {M} }$ . This is not the case of rotating structures, as under the effect of rotation, the damping matrix adopts the form:

(5.26)

Where ${\textstyle {\mathbf {D} _{\,{\hbox{rot}}}}}$ , that accounts for the Coriolis effect, cannot be expressed in terms of ${\textstyle \mathbf {K} }$ nor ${\textstyle \mathbf {M} }$ and neither is a symmetric matrix. Therefore, modal analysis in rotating structures becomes impracticable unless Coriolis effect is neglected. Nonetheless, it is a powerful tool that allows extraction of important information about the system, such as natural frequencies and mode shapes. Moreover, modal decomposition leads to exact results if all the ${\textstyle \mathbf {n} }$ terms of the summation are taken into consideration. This exact solution will help assessing the validity of the approximated results obtained through the numerical methods presented in the next section.

5.2 Numerical integration

In order to overcome the limitations modal analysis is subject to, numerical integration methods arose. The key is to solve the elastodynamic differential equation by integrating over the time domain. As these integrals do not have exact close solution, they are approximated trough a set of numerical methods and algorithms. In particular, the implicit Newmark ${\textstyle {\boldsymbol {\beta }}}$ -method [12] is the one covered in this report. To start with, the time interval is divided into ${\textstyle {\textrm {N}}}$ time steps:

(5.27)

Where ${\textstyle {\textrm {t}}^{n}}$ is a generic time value. Following this description, equation 5.3 is discretised:

(5.28)

Recall that ${\textstyle \mathbf {D} }$ and ${\textstyle \mathbf {K} }$ are no longer constant, as angular acceleration is now taken into consideration. Newmark ${\textstyle \beta }$ -methods are a family of integration schemes in which displacements and velocities at time ${\textstyle {\textrm {t}}^{n+1}}$ are uploaded through the following formulae:

(5.29)

(5.30)

(5.31)

Where ${\textstyle \Delta {\textrm {t}}={\textrm {t}}^{n+1}-{\textrm {t}}^{n}}$ is the time step and ${\textstyle \beta }$ and ${\textstyle \gamma }$ are numerical parameters. The parameter ${\textstyle \gamma }$ introduces damping to the numerical method, controlling what is called artificial viscosity [12]. It is used to reduce noise in the solution. If ${\textstyle \gamma ={\frac {1}{2}}}$ , no damping is added, whereas if ${\textstyle \gamma >{\frac {1}{2}}}$ the artificial viscosity added to the system is proportional to ${\textstyle \gamma -{\frac {1}{2}}}$ . On the other hand, the ${\textstyle \beta }$ parameter defines how much implicit or explicit the method is. Continuing with the formulation of the numerical scheme, if equation 5.31 is introduced into 5.28:

(5.32)

Acceleration at ${\textstyle {\textrm {t}}^{n+1}}$ is obtained from equation 5.30:

(5.33)

Then, expression 5.32 transforms into:

({\tilde {\mathbf {M} }}^{n+1}+\mathbf {K} ^{n+1})\;\mathbf {d} ^{n+1}=\mathbf {F} ^{n+1}-\mathbf {D} ^{n+1}\;{\tilde {\mathbf {v} }}^{n}+{\tilde {\mathbf {M} }}^{n+1}\;{\tilde {\mathbf {d} }}^{n}\;\;,\;\;\;{\hbox{with}}\;\;\;\;\;{\tilde {\mathbf {M} }}^{n+1}={\frac {\mathbf {M} +\gamma \Delta {\textrm {t}}\;\mathbf {D} ^{n+1}}{\beta \,\Delta {\textrm {t}}^{2}}}

(5.34)

The algorithm for the implicit case ( ${\textstyle \beta >0}$ ) can be summarised as:

Compute ${\textstyle {\tilde {\mathbf {d} }}^{n}}$ and ${\textstyle {\tilde {\mathbf {v} }}^{n}}$ through equations 5.30 and 5.31.
Obtain ${\textstyle {\mathbf {d} }^{n+1}}$ solving the system of equations stated in 5.34.
Update accelerations ${\textstyle {\mathbf {a} }^{n+1}}$ and velocities ${\textstyle {\mathbf {v} }^{n+1}}$ using expressions 5.33 and 5.31.

Figure 5.4: Newmark

-method algorithm schematisation for the case of rotating structures.

This scheme is found to be unconditionally stable for ${\textstyle \beta \geq {\frac {\gamma }{2}}\geq {\frac {1}{4}}}$ . It is found that ${\textstyle \beta ={\frac {1}{4}}}$ and ${\textstyle \gamma ={\frac {1}{2}}}$ yield to reasonable convergence rates. If those particular values of numeric parameters are chosen, then the scheme is known as constant average acceleration method. Recall that if ${\textstyle \;\mathbf {F} ^{n+1}=f({\mathbf {d} }^{n+1},\,{\textrm {t}}^{n+1})}$ , then 5.34 would be a nonlinear system of equations, which is to be tackled through numerical schemes such as Newton-Raphson. Fortunately, this is not the case of study and displacements are simply computed as:

(5.35)

5.2.1 Timestep selection

The presented method is a practical tool that allows the solving of complex nonlinear systems. However, convergence of the solution is only granted for sufficient small values of ${\textstyle \Delta {\textrm {t}}}$ . Its value is directly linked with the highest frequency involved in the problem. According to the Nyquist-Shannon theorem, sampling frequency ${\textstyle f_{s}}$ needs to be at least twice the maximum frequency of the signal. In the case of study, two different phenomena are to be taken into account: vibration in the noninertial frame and rotation in the inertial one.

The period of the rotation motion, supposing constant angular acceleration ${\textstyle \alpha }$ and initial velocity ${\textstyle v_{0}}$ , can be posed as:

(5.36)

As rotations undergoing angular acceleration for a long period of time are not usual, let's consider a rotation at very high constant spin speed, say ${\textstyle \Omega =15000\;{\hbox{rpm}}}$ . In consequence, the sampling frequency should be greater than ${\textstyle 500\;Hz}$ , and the timestep smaller than 2 ${\textstyle ms}$ . To put these numbers in context, recall natural frequencies from Figure 5.1: the third mode presented a frequency of ${\textstyle 637\;Hz}$ , implying a sampling frequency greater than ${\textstyle 1275\;Hz}$ .

In general, vibration frequencies are higher than the rotation velocity, and thus the selection of ${\textstyle \Delta {\textrm {t}}}$ will not be influenced by the velocity of the rotation motion. The criteria to follow when picking a sampling frequency is to ensure that it is at least twice the highest vibration frequency to capture.

5.2.2 Validation

In order to asses whether the numerical scheme is performing properly, numerical results will be compared with those obtained with modal analysis. As the latter is restricted to free undamped static structures, the case of study will involve the free vibration of the 2D cantilever rectangular beam presented in 4.7. Imagine the beam under the effect of a distributed load along its upper side, when all of the sudden the force vanishes and the structure starts to oscillate from its deformed configuration around its equilibrium position.

Although highly visual, animations are not a proper way of comparing quantitatively the two responses. Instead, the plot of dynamic displacements at different degrees of freedom will be compared for both modal and numerical cases. In fact, the nodes whose vertical displacement will be graphed are those located at the centre and tip of the mean line of the beam. The modal solution will be obtained taking into account the six first mode shapes with the lowest frequency, whereas the numerical scheme will be based on the constant average acceleration method ( ${\textstyle \gamma ={\frac {1}{2}}}$ and ${\textstyle \beta ={\frac {1}{4}}}$ ). The selected timestep depends of the highest natural frequency to capture, which in this case is linked with the pure bending mode. Hence, ${\textstyle \Delta {\textrm {t}}=1\,ms}$ . The results are compared in Figure 5.5, where dots represent the approximated numerical solution.

Figure 5.5: Comparison of modal and numerical methods to solve the free undamped vibration of a 2D fixed-free beam.

As can be spotted, amplitude in the tip is higher than in the centre. Furthermore, the frequency of vibration is that of the bending motion ( ${\textstyle \sim 41\,Hz}$ ), as was expected from Figure 5.2. The main problem of comparing both solutions is that numerical methods, unlike modal analysis, cannot filter high frequency modes. Thus, the numeric solution is always affected by high frequency components and is not capable of representing the bending motion in isolation. In consequence, modal and numeric results will only be coincident if all the mode shapes are taken into consideration and the time step is sufficiently small. Figure 5.6 compares horizontal displacements at the tip using 600 different eigenvalues, and a numeric timestep of ${\textstyle 0.1\mu s}$ .

Figure 5.6: Modal and numeric (red crosses) methods comparison, this time with higher accuracy.

5.2.3 Limitations

Numerical methods allow complex and even nonlinear systems to be tackled, obtaining results with a proper degree of accuracy. Even so, one of their main drawbacks is linked with the requirement of really small time steps in order to capture structural vibrations. This implies high computation time and memory usage. Recall that in order to animate the results, displacements need to be computed for each timestep. For instance, the displacements of a 3D system with 30000 nodes studied during a timespan of ${\textstyle 10s}$ with ${\textstyle \Delta {\textrm {t}}=0.1\,ms}$ would take up 67 gigabytes of memory. In reality, postprocessing tools also require from other inputs, increasing even more the memory consumption. This makes numeric simulations of complex structures only feasible for institutions with high computational resources.

Another disadvantage of numerical schemes with respect to modal analysis is that no mode shapes nor natural frequencies are obtained. In fact, what is computed through numerical methods is the temporal evolution of the individual displacements at each degree of freedom. The whole phenomena of vibration is difficult to capture by only looking at individual graphs, so postprocessing tools use this information to create an animation of the dynamic response. But animations are not an outright solution, as no information regarding the predominating modes of the structure and its frequency are obtained.

When very different frequencies play an important role in the system, the latter has to be studied at different time scales as there will be motions with very long period coexisting with other ones which elapse very fast. The study of both motions at the same time based on animations is very tricky, and requires from patience and experience. In order to overcome this limitation, the next section reviews a promising method that captures the predominating modes and frequencies of a structure using only the raw data obtained from the numerical approach.

5.3 Singular value decomposition

The singular value decomposition, commonly referred as SVD, is one of the most versatile tools of linear algebra. It allows the factorisation of a rectangular matrix ${\textstyle \mathbf {A} }$ of the form:

(5.37)

where:

${\textstyle \mathbf {A} }$ is a ${\textstyle {\textrm {m}}\times {\textrm {n}}}$ matrix.
${\textstyle \mathbf {U} }$ is a ${\textstyle {\textrm {m}}\times {\textrm {m}}}$ orthonormal matrix known as the matrix of left side vectors (LSV).
${\textstyle \mathbf {V} }$ is a ${\textstyle {\textrm {n}}\times {\textrm {n}}}$ orthonormal matrix known as the matrix of right side vectors (RSV).
${\textstyle {\boldsymbol {\Sigma }}}$ is a ${\textstyle {\textrm {m}}\times {\textrm {n}}}$ diagonal matrix known as the matrix of singular values.

The diagonal entries of ${\textstyle {\boldsymbol {\Sigma }}}$ are referred as singular values ${\textstyle \sigma _{i}}$ , fulfilling that ${\textstyle \sigma _{i}>0}$ . Quite often, the matrix of singular values is sorted in ascending order, so that ${\textstyle \sigma _{i+1}\leq \sigma _{i}}$ . If ${\textstyle \mathbf {A} }$ is a real matrix, then both ${\textstyle \mathbf {U} }$ and ${\textstyle \mathbf {S} }$ are orthonormal: ${\textstyle \mathbf {U} ^{T}\,\mathbf {U} =\mathbf {V} ^{T}\,\mathbf {V} =\mathbf {I} }$ . An intuitive approach to SVD is to think of ${\textstyle \mathbf {A} }$ as a linear transformation. Then, the singular value decomposition can be though as a composition of three geometrical transformations: a rotation ${\textstyle \mathbf {U} }$ , a scaling or stretching ${\textstyle {\boldsymbol {\Sigma }}}$ and another rotation ${\textstyle \mathbf {V} ^{T}}$ .

The SDV can be interpreted as a generalisation of the eigendecomposition, which is restricted to positive semidefinite matrices. In fact, the SVD can be split into a pair of eigenproblems:

{\begin{matrix}\mathbf {A} \;\mathbf {A} ^{T}=\mathbf {U} \,{\boldsymbol {\Sigma }}\,{\overset {\mathbf {I} }{\overbrace {\mathbf {V} ^{T}\,\mathbf {V} } }}\,{\boldsymbol {\Sigma }}^{T}\,\mathbf {U} ^{T}\;\;\;\Rightarrow \;\;\;\mathbf {A} \;\mathbf {A} ^{T}\mathbf {U} =\mathbf {U} \,{\boldsymbol {\Sigma }}^{2}\\\mathbf {A} ^{T}\mathbf {A} =\mathbf {V} {\boldsymbol {\Sigma }}^{T}\,{\overset {\mathbf {I} }{\overbrace {\mathbf {U} ^{T}\,\mathbf {U} } }}\,{\boldsymbol {\Sigma }}\,\mathbf {V} ^{T}\;\;\;\Rightarrow \;\;\;\mathbf {A} ^{T}\mathbf {A} \,\mathbf {V} =\mathbf {V} \,{\boldsymbol {\Sigma }}^{2}\end{matrix}}

(5.38)

Hence, ${\textstyle \mathbf {U} }$ can be thought as the set of orthonormal eigenvectors of ${\textstyle \mathbf {A} \;\mathbf {A} ^{T}}$ and ${\textstyle \mathbf {U} }$ as the set of orthonormal eigenvectors of ${\textstyle \mathbf {A} ^{T}\,\,\mathbf {A} }$ . The singular values of ${\textstyle \mathbf {A} }$ are in fact the square root of the eigenvalues of both ${\textstyle \mathbf {A} \;\mathbf {A} ^{T}}$ and ${\textstyle \mathbf {A} ^{T}\,\,\mathbf {A} }$ . This property brings to light the capacity of the SVD to discover some of the same kind of information as with the modal analysis. In fact, the simplest numerical scheme to tackle the SVD problem consists in solving one of the above eigenproblems, obtaining either ${\textstyle \mathbf {U} }$ or ${\textstyle \mathbf {V} }$ as eigenvectors and ${\textstyle {\boldsymbol {\Sigma }}^{2}}$ as eigenvalues. Then, the remaining matrix of vectors is obtained using equation 5.37.

5.3.1 Truncated SVD

The SVD is an interdisciplinary technique used in diverse fields such as signal processing, pattern recognition, data compression, weather prediction and social-media algorithms among many others. A popular application is related with dimensionality reduction: a process through which a given matrix ${\textstyle \mathbf {A} }$ is approximated by a lower rank matrix ${\textstyle \mathbf {A} ^{(r)}}$ . The SVD provides the best k-rank approximation to ${\textstyle \mathbf {A} }$ of the form:

(5.39)

The key is to perform a SVD picking only the ${\textstyle {\textrm {r}}}$ -largest values of ${\textstyle \sigma _{i}}$ . If ${\textstyle r<n}$ , then the dimension of the matrices involved in the SVD problem are reduced: ${\textstyle \mathbf {U} }$ becomes ${\textstyle {\textrm {m}}\times {\textrm {r}}}$ , ${\textstyle \mathbf {V} ^{T}}$ is now ${\textstyle {\textrm {r}}\times {\textrm {m}}}$ and ${\textstyle {\boldsymbol {\Sigma }}}$ is reduced into a ${\textstyle {\textrm {r}}\times {\textrm {r}}}$ matrix. The truncated version of the SVD is tackled as an optimization of the Frobenius norm of the difference between the original matrix and its approximation. Numerical schemes use a given value of tolerance to limit the truncation error.

Even considering the truncated model, computing the SVD can be extremely time-consuming for large-scale problems. One of the most used techniques to speed up the algorithm relies on block matrix decomposition. The factorisation no longer processes the whole matrix, but performs successive SVDs over smaller sub-matrices. Moreover, in this way the quality of the approximation can be controlled blockwise. These kind of method is known as partitioned truncated SVD.

To improve even more the performance of the algorithm, randomised techniques [10] are used to compute an approximation to the range of ${\textstyle \mathbf {A} }$ . These techniques are of special interest when the data comes from discretization of continuum systems. An intuitive approach is to estimate the range of ${\textstyle \mathbf {A} }$ by picking a random collection of vectors of the matrix and examining the subspace formed by the action of ${\textstyle \mathbf {A} }$ on each of them. For sure, randomised SVD algorithms are not that simple and require a deeper knowledge on linear algebra.

5.3.2 SVD in structural analysis

In the case of study, which is structural analysis, it is from interest to work in the basis defined by the mass matrix ${\textstyle \mathbf {M} }$ instead of the euclidean space. The goal is to re-express ${\textstyle \mathbf {A} }$ in such a way that redundancies are filtered out and actual dominant displacement fluctuation patterns are revealed [13]. In order to develop this normalisation, ${\textstyle \mathbf {M} }$ needs to be a symmetric positive definite matrix such that admits the Cholesky factorisation:

(5.40)

Where ${\textstyle {\overline {\mathbf {M} }}}$ is an upper triangular matrix. Then, the truncated SVD is not performed to the matrix ${\textstyle \mathbf {A} }$ but to the normalised matrix ${\textstyle {\overline {\mathbf {A} }}={\overline {\mathbf {M} }}\;\mathbf {A} }$ :

(5.41)

Finally, the right side vectors are obtained from ${\textstyle {\overline {\mathbf {U} }}}$ :

(5.42)

Following this approach, the SVD can be posed in terms of structural analysis. The matrix ${\textstyle \mathbf {A} }$ is nothing else but the matrix of displacements ${\textstyle \mathbf {d} }$ defined in equation 5.14. Each column contains the displacements of the unconstrained degrees of freedom at a given value of time, whereas each row accounts for the temporal evolution of the displacement at a single degree of freedom. Hence, defining ${\textstyle {\textrm {n}}_{\,\mathbf {l} }}$ as the number of unconstrained degrees of freedom and ${\textstyle {\textrm {n}}_{\,\Delta {\textrm {t}}}}$ as the number of time-steps into which the time domain is split, ${\textstyle \mathbf {A} }$ is a ${\textstyle {\textrm {n}}_{\,\mathbf {l} }\times {\textrm {n}}_{\,\Delta {\textrm {t}}}}$ matrix.

From the structural analysis perspective, ${\textstyle \mathbf {U} }$ is a ${\textstyle {\textrm {n}}_{\,\mathbf {l} }\times {\textrm {r}}}$ matrix containing in each column a certain combination of nodal displacement patterns which are predominating in the system as time passes by. These can be thought, in analogy with modal decomposition, as the principal modes of the structure. In reality, left side vectors are not modal mode shapes but a certain combination of the sadistically reigning ones. Hence, on now on the vectors ${\textstyle \mathbf {U} _{i}}$ will be referred as singular modes.

On the other hand, ${\textstyle \mathbf {V} ^{T}}$ is a ${\textstyle {\textrm {r}}\times {\textrm {n}}_{\,\Delta {\textrm {t}}}}$ matrix each row of which contains a statistically predominant temporal evolution of the whole structure. In contrast with modal analysis, each ${\textstyle \mathbf {V} _{i}}$ describes a vibration motion which is not necessarily harmonic. The vectors ${\textstyle \mathbf {V} _{i}}$ will referred as singular oscillations.

The grace of this method is that each ${\textstyle \mathbf {V} _{i}}$ is linked with ${\textstyle \mathbf {U} _{i}}$ , allowing a certain intuition on how each of the singular modes evolves in time. Through an adequate analysis, one can find out the predominating frequencies of each singular oscillation (RSV), and consequently the vibration frequencies ${\textstyle {\overline {\boldsymbol {f}}}_{i}}$ associated to each singular mode (LSV). Moreover, singular values ${\textstyle \sigma _{i}}$ give a hint on the intensity of each singular mode, allowing the truncated study of the system accounting only for the most relevant ones. The individual motions can then be animated using GID, obtaining a more clear sight of the phenomena than the one obtained with numerical integration.

As it has been pointed out, singular oscillations are no longer a pure harmonic and thus may contain the contribution of more than one natural frequency. In order to extract information regarding predominant frequencies, two numerical tools can be used. The most intuitive approach is based on looking for the local maxima of ${\textstyle \mathbf {V} _{i}}$ , and approximating the natural frequency as the average of the time intervals between successive maximums. The main drawback of the latter method is that if more than two frequencies are involved in the temporal description of the motion, only the highest one can be extracted from ${\textstyle \mathbf {V} _{i}}$ . Furthermore, the induced error on the measure of the latter is no longer negligible.

A more feasible approach consists in applying the discrete Fourier transform on ${\textstyle \mathbf {V} _{i}}$ using in-built ${\textstyle MATLAB}$ functions. By doing so, the singular oscillation is mapped into the frequency domain. Predominating frequencies of each singular mode can be found by looking for the peaks on the power spectral density estimation constructed from the Fourier transform. Recall that the obtained vectors ${\textstyle \mathbf {U} _{i}}$ and ${\textstyle \mathbf {V} _{i}}$ have been scaled — in space and temporal coordinates, respectively — during the SVD, and thus information regarding the sampling frequency used during the numeric integration is needed to recover the actual frequencies involved in the vibration motion.

For sure, the development of a numerical algorithm capable of handling the truncated SVD problem in a short amount of time is far out of the scope of this report. Instead, a ${\textstyle MATLAB}$ tool developed by prof. Hernández J. will be used to perform a partitioned randomised SVD. This software limits the number of terms of the truncated approximation ${\textstyle \mathbf {A} ^{(r)}}$ by using a value of tolerance introduced by the user.

5.3.3 Validation

In order to determine whether the singular value decomposition is a suitable technique for structural analysis, its performance will be assessed by comparing the singular modes and predominant frequencies with the results obtained through modal analysis. Again, the case of study involves a 2D cantilever beam in a free undamped oscillation, as covered in 5.5. Let's consider the dynamic displacements obtained from modal decomposition as the input data matrix for the SVD. The goal is to check whether this method is able to recover the information regarding natural modes and frequencies properly.

In this particular case, the modal response has been approximated using only the six natural modes with the lowest frequency. The sampling frequency used to discretise equation 5.13 is ${\textstyle 500000\,Hz}$ . On the other hand, the truncated SVD will only account for the six more relevant singular modes. To help comparison with modal analysis, one may define the normalised singular values as

(5.43)

With this in mind, the SVD and modal results are compared in the following table:

Table 5.1. Frequencies and intensities comparison: SVD vs. modal analysis

The above table reveals an outstanding degree of accuracy in the frequency determination using the SVD. The singular modes associated to the four lower predominant vibration frequencies are displayed in Figure 5.7.

Figure 5.7: Singular modes associated to the first four lower vibration frequencies.

Comparing the above illustration with natural modes (Figure 5.1), one finds out that they are recovered almost exactly. The only exception is found in the third natural mode, related with the axial deformation. In this case, neither the predominant frequency nor the singular mode match exactly with the modal results. The reason behind this discordance is that the axial mode is merely activated in a bending scenario. In consequence, being the SVD a statistical method, less reliability is to be expected when recovering the axial mode.

With all, the singular value decomposition has been proven to be a good approach in recovering relevant information from a given matrix of data. In structural analysis, it is a promising tool for obtaining the predominant motions when modal analysis is not practicable. The next chapter will cover the application of the SVD to a set of rotating scenarios.

6. Analysis of rotating beams using the SVD

Once the methodology for handling rotating structures in non-stationary regime is fully developed, it is time to tackle a few scenarios involving rotating machinery. The simulations will be carried taking advantage of the Newman- ${\textstyle \beta }$ algorithm presented in the previous section, and then the SVD will used to disclose the dominant modes and frequencies of the dynamic response. The results will be compared with those obtained through modal analysis, which is restricted to non-rotating vibrations. Hence, a hint on the effects centrifugal, tangential and Coriolis forces have on rotating structures will be brought to light in this chapter.

Postprocessing of results will be carried out by GiD, which uses the generic spatial basis ${\textstyle \left\{{\textrm {xyz}}\right\}^{T}}$ . The simulations performed in the present chapter will consider only variables in the local frame of reference ${\textstyle \left\{\varsigma \,\varrho \,\varpi \right\}^{T}}$ , so although some figures incorporate a set of ${\textstyle {\textrm {xyz}}}$ axes, keep in mind that they represent in fact the corotational frame.

6.1 Two-dimensional cantilever beam

The first case of study continues with the analysis of the rectangular 2D cantilever beam studied in the previous section. The beam is ${\textstyle 0.2\times 2\,m}$ and modelled using 1000 quadrilateral elements. The main difference is that the beam now undergoes a rotation perpendicular to the plane, with its axis located at the centre of the fixed end of the beam. The beam is made from 6061 T6 aluminium alloy ( ${\textstyle E=70\,GPa}$ , ${\textstyle \nu =0.3}$ ), and the considered density is ${\textstyle \rho =\;2700kg/m^{2}}$ .

Figure 6.1: Rectangular cantilever beam subject to rotation.

The loading scenario is quite similar as the presented in 4.6: the beam is under the effect of a distributed ${\textstyle \varrho }$ load on its upper side (-20 kN/m), while it is rotating at a certain angular velocity. Initial displacements under these conditions will be computed using the static FEM software. The dynamic case will account for a sudden vanishing of the distributed load and the consequent vibration of the beam. This oscillating behaviour will be analysed under the effect of different angular speeds and accelerations. No internal damping will be considered in this section. For comparison purposes, the first eight natural frequencies and modes of the presented structure can be found in the report attachments (Figure D.1.1).

6.1.1 Constant angular velocity

The first analysis will involve constant angular speed. The goal is to keep track of the vibration frequencies at different angular velocities, and compare the results with the softening effect predicted by modal analysis. The imposed angular speed cannot be indiscriminately high, as deformations would no longer be considered small. Using the 1D solver developed in 4.2.2, it has been found that elastic displacements reach a 1% of the length of the beam when ${\textstyle \Omega =4000}$ rpm. Hence, this value will not be exceeded in this analysis.

It is important to select an adequate timespan and timestep, as numerical results are sensitive on these parameters. The timestep is chosen depending on the maximum frequency to capture, in this case ${\textstyle f_{\hbox{max}}\approx 2500\,Hz}$ , corresponding to the eighth natural frequency. Hence, ${\textstyle \Delta \,{\textrm {t}}<0.2\,ms}$ . In order to gain accuracy, the used timestep will be ${\textstyle \Delta \,{\textrm {t}}=5\,\mu s}$ . To keep computation time low, the timespan will be no longer than ${\textstyle 0.2\,s}$ . With these numeric parameters, and using the constant average acceleration method, the mean time of computation is around 25 minutes. Smaller time-steps or greater time-spans would be too much time-consuming.

Campbell diagram

The dynamic analysis is to be developed for three different angular velocities: 50, 100 and 200 rad/s. The SVD is truncated to the eight higher singular values. The left and right side vectors obtained for each case are displayed in D.1.2. To ease visualisation of the RSVs, each one is plotted for a certain fraction of the total timespan. The dominant frequencies are found looking for the peaks of the discrete Fourier transform of each RSV. The obtained frequencies are presented in Figure 6.2, and compared with the softening effect predicted by modal decomposition.

Figure 6.2: Campbell diagram: comparison between modal prediction and SVD results.

As can be spotted, in most cases the dominant frequency obtained using SVD matches almost exactly with the predicted natural frequency. In other cases (modes 5, 6 and 8), the values are not coincident but very similar ( ${\textstyle \pm 1Hz}$ ). In all likelihood, these differences can be attributed to the probabilistic behaviour of the SVD, as the numeric sample has a timespan of only ${\textstyle 0.2\,s}$ . Even so, in all cases the vibration frequencies predicted by the SVD decrease with the spin speed, matching with the expected softening effect. Thereby, it can be deduced that neither Coriolis nor centrifugal forces induce a significant change in the natural frequencies of rotating beams. In the same way, by comparing the LSVs from D.1.2 with the natural modes of the beam, it can be stated that mode shapes do not either change with the angular velocity.

Recall that the obtained results shown in D.1.2 are just predominant patterns of the dynamic response. Thus, the RSVs representing reigning vibrations are no longer pure harmonics, but mostly the result of a certain combination of natural frequencies. In most cases, they match almost perfectly to a sine or cosine function, and hence the predominant frequency is easy to determine. That is not the case of the RSVs associated to axial vibrations (natural modes 2 and 7 from Figure D.1.1), which are the result of more than one predominant frequency. For instance, the pure axial motion has a ${\textstyle 0}$ Hz frequency associated to it, caused by the offset in displacements the centrifugal force induces. To visualise the frequencies involved in the vibration phenomena, Figure D.11 displays the power spectral density estimation of the sum of the right side vectors.

Singular values

Although Coriolis and centrifugal acceleration induce little changes in natural frequencies, they do change considerably the intensity in which natural modes show up. Considering the normalised singular values as the indicator of intensity, Figure D.8 shows how axial modes (2 and 7) gain intensity as rotation speed increases, whereas flexural modes weaken. This effect is caused by the centrifugal force, the responsible of axial excitation, which increases proportionally to ${\textstyle \Omega ^{2}}$ . The numeric values of predominant frequencies and normalised singular values are presented in table D.2.

Reactions

Another interesting variable to keep track of is the reaction torque, which at last instance is the one determining the required torque to spin at the desired speed. For the presented scenario and studied angular speeds, the reaction torque is found to oscillate along with the beam, as displayed in Figure D.9.

As the motion is not under the effect of angular acceleration, no rigid body reactions are present. Thus, Figure D.9 only accounts for the elastic reactions described in 3.21. The torque is mainly induced by the effect of transverse reactions, as local horizontal reactions due to centrifugal forces produce no net torque around the specified rotation axis. Hence, the required torque is closely related with the bending modes, and thus, with its natural frequency, which decreases with ${\textstyle \Omega }$ . In this case, reactions are not generated by a certain applied force but as a consequence of the inertia forces induced by the high accelerations present in the vibration phenomena.

At some intervals, the required torque becomes negative: half of the time, the beam is bending in the direction of rotation, tending to accelerate the structure. Hence, the overall work is zero. This is an artificial scenario: the angular speed of actual rotors is never constant, whereas the torque injected by the driving device is in fact uniform. In real cases, the vibration of the beam accelerates and decelerates the rotation motion in a small amount (as stated in 1.2). Moreover, real rotors work within viscous fluids, which act as dampers of the vibration motion.

Displacements

Regarding dynamic results, displacements at the tip of the beam will be reviewed for different rotation velocities. Displacements in the local perpendicular direction are induced mostly by the pure flexural mode, and hence its vibration motion is easy to capture as it is associated to a low frequency. The amplitude of these displacements does not increase significantly with ${\textstyle \Omega }$ , but it is the decrease in frequency the most remarkable effect rotation velocity introduces. Figure D.10 presents the numeric results and their lower rank approximation computed using the truncated SVD with eight singular values.

On the other hand, axial displacements are highly dependant on ${\textstyle \Omega }$ , as the centrifugal force that excites axial modes increases with the square of the angular velocity. In fact, the amplitude increases by a factor of 16 when comparing Figure D.12 (a) with Figure D.12 (c). Moreover, the local horizontal displacements are linked to a relatively high frequency, and thus the truncated SVD is not capable to recreate the numerical result exactly. This issue is easily overcome by decreasing the value of tolerance of the SVD.

Stresses

In structural analysis, one of the most usual parameters used to determine whether the structure is going to fail or not is the security factor. In this example, computation of stresses is quite straightforward, as no internal damping is considered. Stresses and strains are obtained directly from the numerical scheme, but unlike displacements, they are not passed into the SVD for dominant patterns recognition. When dealing with stresses, their distribution is not so important as it is the peak value, from which depends the security factor. Figure 6.3 displays the peak values of ${\textstyle \sigma _{\varsigma \varsigma }}$ , ${\textstyle \sigma _{\varrho \varrho }}$ and ${\textstyle \tau _{\varsigma \varrho }}$ as time passes by for the most favourable case ( ${\textstyle \Omega =50}$ rad/s).

From Figure 6.3, it is found that maximum stresses oscillate at twice the natural frequency of the flexural mode. That is caused by the fact that maximum stresses take place at maximum deformed configurations — highest strains — which occur two times per oscillation. For the used aluminium, the yield strength is around 276 MPa [14], and so the security factors are 50, 13, and 4 for 0, 50 and 100 rad/s, respectively.

Figure 6.3: Temporal evolution of the maximum stresses (absolute value) for

= 50 rad/s.

Results animation

To close down this section, the outputs from the dynamic analysis will be animated through ${\textstyle \mathbf {GiD} }$ . However, postprocessing is a heavy task for computers with humble computation capabilities, as results are encoded in ASCII. For that reason, only the case corresponding to ${\textstyle \Omega =50}$ rad/s will be dynamically animated. To simplify the animation, the deformed configuration will be displayed in the rotating frame, and only stresses in the local horizontal direction will be represented through a colour map.

During the bending motion, one expects to found a symmetric distribution of stresses with respect to the beam's axis: one side under traction and the other under compression. Figure D.13 shows how this symmetry is broken in rotating scenarios, as the centrifugal force introduces an axial displacement proportional to the distance from the rotation axis that counteracts the compressing stresses.

6.1.2 Constant angular acceleration

The second and last analysis regarding the rectangular cantilever beam accounts for the effect of a constant acceleration. Recall that angular acceleration was neglected in modal analysis, as variable rotation speed is not compatible with its formulation. Hence, this analysis aims to unfold the effect angular acceleration has on rotating structures. For the sake of simplicity, only two simulations will be carried out, corresponding to ${\textstyle \alpha =5}$ and ${\textstyle 50\,{\hbox{rad/s}}^{2}}$ . In both cases, the initial angular velocity will be ${\textstyle \Omega _{0}=50\;{\hbox{rad/s}}}$ .

The numerical parameters will be the exact same as in the case of constant spin speed. Thus, the timespan will not be large enough for the angular acceleration to introduce significant changes in angular velocity — in the case of higher acceleration, ${\textstyle \Delta \Omega =10\;{\hbox{rad/s}}}$ . For this reason, no Campbell diagram will be represented, as the softening effect induced by ${\textstyle \alpha }$ will be relatively small and difficult to measure. Remember that the rotation stiffness matrix is proportional to:

(6.1)

Thus, as ${\textstyle \Omega ^{2}\gg \alpha }$ , no significant effect is expected. The first difference with the previous case is the tangential force introduced by ${\textstyle \alpha }$ , acting perpendicular to the beam. In consequence, the flexural mode is highly activated when angular acceleration is present. In fact, the flexural motion is so predominant that the SVD is no longer capable of extracting neither the natural frequency nor the mode of the axial vibration. Figure D.16 reveals how the second predominant mode from SVD now corresponds to a combination of both axial and perpendicular displacements. Intensity of axial modes (three and seven from Figure D.1.1) decreases with angular acceleration, as shown in Figure D.18, that compares normalised singular values of different natural modes for different values of ${\textstyle \alpha }$ .

Notice that the truncated SVD is not capable of recovering the seventh natural mode when ${\textstyle \alpha =50\;{\hbox{rad/s}}^{2}}$ . In fact, in this case the seventh and eight LSVs are the same and correspond to the eight natural mode ${\textstyle \Phi _{8}}$ . Predominant frequencies and singular values obtained are summarised in Table D.3.

Displacements

Again, displacements at the tip will be considered representative of the vibration phenomena. Displacements in the perpendicular direction oscillate at the first natural frequency, but this time the vibration is not around zero but to a new equilibrium position which depends on the magnitude of the applied tangential force. The higher the angular acceleration, the more deformed the configuration at equilibrium will be. This behaviour can be spotted in Figure D.19.

On the other hand, axial displacements are relatively smaller (one order of magnitude below perpendicular ones). They oscillate at two different frequencies, corresponding to the flexural and axial natural modes. It is worth pointing that, as the angular velocity increases with time, so does the centrifugal force. For that reason, the mean value of local horizontal displacements increases with time, as captured by the second RSV in Figure D.17. The rate of this rise is proportional to the angular acceleration. This effect can be identified in Figure 6.4, where horizontal displacements are represented. The increase in length of the beam ${\textstyle \Delta \varsigma }$ during the analysed timespan are ${\textstyle 0.11}$ and ${\textstyle 0.01}$ mm for ${\textstyle \alpha =50}$ and ${\textstyle 5\;{\hbox{rad/s}}^{2}}$ , respectively. If the angular acceleration (for instance ${\textstyle \alpha =50\;{\hbox{rad/s}}^{2}}$ ) persisted in time, that would result in an increase of half millimetre in the beam's length every second.

Figure 6.4: Local horizontal displacements (m) at the tip of the beam for different values of angular acceleration.

Reactions

As in the previous case, the obtained reactions are those the driving device must inject into the structure to ensure the desired acceleration. In this case of constant ${\textstyle \alpha }$ , the elastic reaction torque oscillates at the bending frequency, but with a mean value different from zero as the equilibrium position of the beam is no longer the horizontal configuration. Figure D.20 displays how the mean value of elastic reaction torque increases with ${\textstyle \alpha }$ .

It is worth pointing that, in the case of variable angular speed, a rigid body motion must be taken into account and added to the previous one, as stated in equation 3.2.3. In this case, the rigid body torque is constant:

(6.2)

This constant value is to be added to the elastic torque from D.20, which is significantly smaller. Recall that this is not a real scenario, as the kinematics of the motion are often the unknown rather than the input parameter. The required power the driving engine must provide is given by the product of the total reaction torque (rigid body and elastic) by the instantaneous angular velocity ${\textstyle \Omega }$ .

Stresses

In order to validate structural reliability, one needs to focus on the study of stresses, which determine the security factor. Again, as no internal damping is considered, computation of stresses is done as in the static scenario. For simplicity, only peaks values will be considered in this section. As expected, stresses increase with ${\textstyle \alpha }$ , although they are most sensitive to a change in angular velocity. Figure 6.5 reveals how stresses in the local horizontal direction are the higher ones, and oscillate in the flexural natural frequency. The peaks in stresses correspond to the maximum deformed configurations.

Figure 6.5: Temporal evolution of the maximum stresses (absolute value) for

= 50 rad/s and

= 50 rad/s

.

It is worth pointing that the mean value of ${\textstyle \sigma _{\varsigma \varsigma }}$ increases with time, consequence of the increasing centrifugal force that is acting on the structure. Although only absolute values are represented in Figure 6.5, looking into the dynamic results (Figure D.21) one finds out that traction stresses are the most relevant contribution. As seen in the previous case, centrifugal force tends to pull the beam decreasing the effect of compression stresses.

Results animation

To ease visualisation of results, dynamic displacements and stress distribution will be animated in the rotation frame trough GiD. Only the case of higher angular acceleration ( ${\textstyle \alpha =50}$ rad/s ${\textstyle ^{2}}$ ) is presented in Figure D.21. The animation captures one oscillation of the beam, with a timespan of 25 ms. As commented, higher stresses are found in the most deformed configurations of the beam.

6.2 Three-dimensional cantilever beam

For the next case of study, a new geometry will be defined. In order to keep the problem simple, a cantilever beam will again be the structure to study. In this case, however, the beam is modelled in three dimensions using 3200 hexahedral elements. The beam is ${\textstyle 2\times 0.25\times 0.25\,m}$ , with a thickness of ${\textstyle 0.05\,m}$ . Again, this study considers that the structure undergoes a rotation in the vertical - ${\textstyle z}$ - axis. The material from which the beam is made of is aluminium, with a density ${\textstyle \rho =2700\;Kg/m^{3}}$ , a Young's modulus ${\textstyle E=70\;GPa}$ and a Poisson's coefficient ${\textstyle \nu =0.3}$ .

Figure 6.6: 3D cantilever beam subject to rotation around the vertical axis.

The loading scenario is the following: the beam is rotating at a given angular speed, subjected to ${\textstyle \alpha =10}$ rad/s ${\textstyle ^{2}}$ , when the angular acceleration vanishes instantaneously and so does the tangential force. This section will analyse the vibrations this sudden stop in acceleration induces, for different values of angular speed ${\textstyle \Omega }$ . Again, no source of damping will be considered. In contrast with the last analysis, this time the weight of the beam will be taken into consideration, and supposed constant and negative in the ${\textstyle z}$ axis. For comparison purposes, the first eight natural frequencies and modes of the presented structure can be found in D.2.1.

6.2.1 Constant angular velocity

For the sake of simplicity, the presented geometry will only be studied for the case of constant angular velocity. In fact, the dynamic behaviour will be analysed for three values of ${\textstyle \Omega }$ : 25, 50 and 100 rad/s. The goal is to keep track of the vibration frequencies using the truncated SVD as ${\textstyle \Omega }$ increases. For the numerical computations, the used timestep will be ${\textstyle \Delta \,{\textrm {t}}=25\,\mu s}$ . To keep computation time low, the timespan has been selected in order to capture a full oscillation of first natural mode, so ${\textstyle {\textrm {T}}_{f}=17\,ms}$ . For this case, the truncated SVD will compute the eight highest singular values and their associated LSVs and RSVs.

Campbell diagram

One of the main applications of the truncated SVD is to recover predominant frequencies of the dynamic vibration motion. The goal is to recover these frequencies for different values of ${\textstyle \Omega }$ applying the discrete Fourier transform to each RSV. Recall that in the 2D case, predominant frequencies were compared with the softening effect predicted by modal analysis (see equation 5.25 ). However, this approach is not possible in three dimensions as the rotation stiffness matrix is no longer proportional to the mass matrix. The eigenvalue problem admits no further simplification than:

(6.3)

The only way to keep track of the natural frequencies as angular velocity increases is by solving the previous equation for different values of ${\textstyle {\mathbf {K} _{\hbox{rot}}}\,(\Omega )}$ , obtaining numerically the desired Campbell diagram ( ${\textstyle \omega \;vs.\;\Omega }$ ). By doing so, it is found that in most cases natural frequencies decrease in the same way as they did in 2D. The biggest difference is that in three dimensions, natural frequencies linked with modes acting in the ${\textstyle z}$ direction — the rotation axis — decrease much slower and in a different pattern.

This effect is displayed in a Campbell digram in D.29, where modal prediction and SVD results are compared. It is worth pointing that, in still conditions, natural modes and frequencies appear in pairs, as the beam is symmetric with respect to the ${\textstyle \varsigma \varrho }$ and ${\textstyle \varsigma z}$ planes. However, as rotation velocity increases, vibration frequencies of the modes acting on ${\textstyle \varrho }$ and ${\textstyle z}$ diverge, as their decrease rate is not the same (see ${\textstyle \Phi _{7}}$ vs ${\textstyle \Phi _{8}}$ in D.29).

Recall that one of the advantages of the SVD with respect to modal analysis is that it accounts for the effect of Coriolis, centrifugal and tangential forces. On the other hand, its main drawback is that as it is a probabilistic method, great timespans and small timesteps are needed to get a reliable value of predominant frequency. Nonetheless, the results displayed in D.29 agree ( ${\textstyle \pm 1Hz}$ ) with the predicted softening effect. It is found that the more intensively a mode is activated, the greater the accuracy in which the SVD can extract its predominant frequency. The RSVs from which the predominant frequency has been obtained are plotted for different values of ${\textstyle \Omega }$ in D.24, D.26 and D.28.

Another limitation of the SVD is that it is not able to isolate natural modes. The LSVs the SVD computes are just predominant patterns, which quite often correspond to the combination of two or more natural modes. Figures D.23, D.25 and D.27 — which represent the predominant modes of the beam for different values of ${\textstyle \Omega }$ — evince for instance how the first predominant mode ${\textstyle \mathbf {U} _{1}}$ is in reality the combination of two natural modes: axial in ${\textstyle \varsigma }$ ( ${\textstyle {\boldsymbol {\Phi }}_{6}}$ ) and bending in ${\textstyle z}$ ( ${\textstyle {\boldsymbol {\Phi }}_{2}}$ ). The higher ${\textstyle \Omega }$ is, the higher contribution the axial mode has, as centrifugal forces increase with ${\textstyle \Omega ^{2}}$ whereas weight is constant.

Displacements

Once more, displacements at the tip — upper side, at half the thickness — will be considered representative of the vibration phenomena. This time there is not a predominant direction in which displacements are considerably greater than in the others. Displacements in the ${\textstyle \varrho }$ axis are those with the lowest frequency and linked with the natural mode ${\textstyle {\boldsymbol {\Phi }}_{1}}$ (see D.2.1), as the beam was under the effect of a transverse force which vanished at ${\textstyle t=0}$ , oscillating thus around the equilibrium position ${\textstyle \varrho =0}$ (Figure D.30).

The structure undergoes as well a bending in the vertical direction ${\textstyle z}$ , but with a much higher frequency and a considerably smaller amplitude of oscillation. This effect is produced because, unlike in the ${\textstyle \varrho }$ direction, the distributed force in ${\textstyle z}$ — the weight — does not vanishes but remains constant as time passes by. Hence, vibrations in the vertical axis are not caused by the bending mode ${\textstyle {\boldsymbol {\Phi }}_{2}}$ but are the result of higher frequency modes. Even though vibration amplitudes are relatively small in ${\textstyle z}$ , displacements on this axis have a considerable offset induced by the weight, as can be spotted in Figure D.31.

A similar offset is found in the axial displacements, induced by the constant centrifugal force. As in the vertical direction, displacements in ${\textstyle \varsigma }$ are subject to high frequency vibration, but to small oscillation amplitudes. One of the vibration frequencies corresponds to ${\textstyle \omega _{1}}$ , as bending in ${\textstyle \varrho }$ also causes a small displacement in ${\textstyle \varsigma }$ . For sure, there is as well a component linked with the axial mode ${\textstyle {\boldsymbol {\Phi }}_{6}}$ , but its amplitude is considerably smaller. These three effects — offset, ${\textstyle \omega _{1}}$ and ${\textstyle \omega _{6}}$ — can be easily identified in Figure 6.7.

Figure 6.7: Temporal evolution of

displacements at the tip of a cantilever beam under constant angular velocity.

With no surprise, all displacements increase (in absolute value) with ${\textstyle \Omega }$ , being the ones in ${\textstyle \varsigma }$ the more susceptible to changes in angular velocity, as the centrifugal force that acts in the axial direction is proportional to ${\textstyle \Omega ^{2}}$ .

Reactions

From a dynamic point of view, it is always interesting to keep track of the elastic reactions. Remember that like in the previous case, the beam is considered to rotate at a given angular speed, set by the driving device (recall the assumptions made in 2.1). This implies that no matter how great the deformations are, the device will always inject the necessary torque to ensure the desired kinematic conditions. Although it is an artificial scenario, it allows uncoupling of rigid body and elastic equations, simplifying considerably the dynamic study.

In this particular case, the elastic reaction torque in ${\textstyle z}$ oscillates with the frequency linked with the transverse motion, the most predominant one. The latter result is quite intuitive: the more deformed the configuration is, the more torque is required. When the deformation is in the direction of rotation, the reaction torque becomes negative. This oscillation in the torque is not caused by the effect of any real force, but because of the inertia induced by accelerations. On the other hand, the biggest component of the torque is actually constant, and produced by the centrifugal force. Thereby, the mean value of the elastic reaction torque in the rotation axis increases with ${\textstyle \Omega }$ , as graphed in Figure D.32.

Stresses

Once more, validation of the structure relies in the study of stresses. In order to simplify the analysis, only the maximum absolute value of stress has been dynamically tracked. Without considering any source of internal damping, these values are found to oscillate at twice the frequency of the flexural motion. Recall that stresses and strains adopt its higher value at maximum deformed configuration, which takes place twice per period. With all, it has been found that like in the 2D case, the higher stresses are those in the direction of the centrifugal force — that is, ${\textstyle \sigma _{\varsigma \varsigma }}$ . Maximum stresses for ${\textstyle \Omega =50}$ rad/s are represented in Figure 6.8.

Figure 6.8: Temporal evolution of the maximum stresses (absolute value) for

= 50 rad/s.

In 3D cases where six different components of stress play a significant role, its worth considering the equivalent stress. Using the Von Mises yield criterion for beams and considering a 6061 T6 alloy, it is found that the security factor of the structure is ${\textstyle 14.5}$ .

Results animation

The dynamic results have been postprocessed and animated using GiD. The deformed configuration and the stress map for different timesteps are presented in Figure D.33. For the sake of simplicity, only ${\textstyle \sigma _{\varsigma \varsigma }}$ stresses are coloured. Regarding displacements, it is relatively difficult to keep track of the bending motion in ${\textstyle \varrho }$ , as axial displacements are considerably high and eclipse the flexural vibration.

7. Study of an aeronautics-related case

So far, the presented chapters encompassed a relatively abstract subject, and has only been put in practice to test the capabilities of the SVD in beam analysis. The proposed cases where artificial scenarios that could be easily studied by following the hypotheses posed in 2.1. The aim of this chapter is to modify the developed numerical tool in order to tackle a case of aeronautical interest, which is the starting of an helicopter rotor.

7.1 Case of study

The main goal of the presented chapter is the study of the dynamic behaviour of an helicopter's rotor during the starting of the engines. That is, the blades will start from idle configuration and will end rotating at its design angular velocity. Helicopter dynamics is a wide and complex subject combining both structural analysis and aerodynamics. As the understanding of the whole phenomena is not the goal of the present report, the following hypotheses will be assumed:

Blades will be modelled as two-dimensional structures. Three-dimensional effects (structural and aerodynamic) will be neglected, as well as gravity forces.
Rotation axis is fixed and perpendicular to the plane of the blade.
The fluid will not be simulated but rather aerodynamic forces will be modelled.
The rotor will produce zero lift, and thus, no induced velocity along the span of the beam will be generated. Parasite drag is the only aerodynamic force acting on the rotor.
The helicopter is performing hovering flight with no wind speed.
The blade is fixed-free: no articulated parts are considered. Hence, the mechanical twist angle will always be zero.

The rotor is made of three equally distributed blades, whose geometry is defined in Figure 7.1. The radius of the blade is ${\textstyle R=5\,m}$ , whereas the chord is ${\textstyle 0.2}$ and ${\textstyle 0.4\,m}$ at the root and the tip, respectively. The maximum chord is ${\textstyle 0.55\,m}$ and its located at the 8% of the span. The total surface of the blade is ${\textstyle A_{\,{\hbox{blade}}}=2.335\,m^{2}}$ . The blade has been modelled using a CAD software (SolidWorks) and then exported to GiD to proceed with the mesh generation. Its geometry is defined in Figure 7.1.

Figure 7.1: Geometry and of a 2D blade.

In order to mesh the above geometry using only quadrilateral elements, the blade needs to be modified. By eliminating curves and roundings, the blade is discretised using a structured mesh made of 1340 elements (Figure 7.2).

Figure 7.2: Structured mesh of a simplified 2D blade.

In three dimensions, the blades can be thought as an extrusion of a NACA 0012 airfoil, whose cross section ${\textstyle A_{\,{\hbox{airfoil}}}=0.02081\,m^{2}}$ has been computed integrating the profile shape [15]. To approximate the superficial density of the blade, it will be considered that its total mass is one fifth of the mass of a solid aluminium blade:

(7.1)

Regarding the driving device, it will be supposed that the helicopter has an alternative engine able to inject power to the rotor. As engines cannot produce power when the angular velocity is zero, an electric motor will be used to generate a torque that starts the rotation motion. Remember that the power ${\textstyle {\textrm {P}}}$ of a rotating object is given by:

(7.2)

Where ${\textstyle {\textrm {T}}}$ is the torque produced by the device. To avoid sudden jumps in the variables, neither the torque nor the power will be injected instantaneously. Instead, they will vary from an initial to a final value describing a sinusoidal shape:

{\begin{matrix}{\textrm {T}}={\textrm {T}}_{0}+{\textrm {T}}_{\hbox{nom}}\,\sin \left({\frac {\pi }{2\,{\textrm {t}}_{a}}}\,{\textrm {t}}\right)\qquad 0\leq t<{\textrm {t}}_{a}\,,\qquad \qquad {\textrm {T}}={\textrm {T}}_{\hbox{nom}}\qquad {\textrm {t}}_{a}\leq t<{\textrm {t}}_{b}\\\\{\textrm {P}}={\frac {1}{2}}\left({\textrm {P}}_{0}+{\textrm {P}}_{\hbox{nom}}\right)+{\frac {1}{2}}\left({\textrm {P}}_{\hbox{nom}}-{\textrm {P}}_{0}\right)\,\sin \left({\frac {\pi }{{\textrm {t}}_{c}-{\textrm {t}}_{b}}}\,{\textrm {t}}+\pi \left({\frac {1}{2}}-{\frac {{\textrm {t}}_{c}}{{\textrm {t}}_{c}-{\textrm {t}}_{b}}}\right)\right)\qquad {\textrm {t}}_{b}\leq t<{\textrm {t}}_{c}\\\\{\textrm {P}}={\textrm {P}}_{\hbox{nom}}\qquad \qquad {\textrm {t}}_{c}\leq t\end{matrix}}

(7.3)

${\textstyle {\textrm {W}}_{0}}$ and ${\textstyle {\textrm {T}}_{\hbox{nom}}}$ are the initial and nominal values, respectively, of torque generated by the electric motor. The motor reaches the nominal value of torque at ${\textstyle {\textrm {t}}_{a}}$ , and maintains it constant until ${\textstyle {\textrm {t}}_{b}}$ , when the engine starts operating. To avoid discontinuities in the injected power, the following statement has to be met:

(7.4)

The engine starts to operate at ${\textstyle {\textrm {t}}_{b}}$ and reaches its design value ${\textstyle {\textrm {P}}_{\hbox{nom}}}$ at ${\textstyle {\textrm {t}}_{c}}$ , maintaining this value constant for further values of time until normal flight conditions are reached (Figure 7.3). No automatic control that regulates the power as a function of the angular velocity is taken into consideration.

7.2 Rigid-elastic coupling

Unlike in previous chapters, kinematics of the rotating motion are no longer known. This implies breaking the first hypothesis defined in 2.1, and thus, the developed model is no longer applicable. In previous cases (see 6), angular velocity and acceleration were known at every time, being the reaction forces the unknowns. In this case, what is known is the power provided by the engine, whereas the unknowns are the angular variables ${\textstyle \theta }$ , ${\textstyle \Omega }$ and ${\textstyle \alpha }$ . This is a major change in the problem formulation, as rotation and vibration motions become coupled. From solid mechanics, it is found that the relationship between the torque and the angular velocity of a rigid body is given by:

(7.5)

Where ${\textstyle a}$ is a generic point where the torque is to be evaluated, ${\textstyle \mathbf {J} _{a}}$ is the inertia tensor of the solid with respect to ${\textstyle a}$ , ${\textstyle {\boldsymbol {\Omega }}}$ is the vector of angular velocities, ${\textstyle \mathbf {a} _{a}}$ the linear acceleration of the point ${\textstyle a}$ and ${\textstyle \mathbf {AG} }$ the vector pointing from ${\textstyle a}$ to the centre of mass of the solid. If the inertia matrix is posed in such a way that is constant in time, the above expression transforms into Euler's equations for rigid bodies (see 1.2). Luckily enough, with the assumed hypotheses, the torque equation can be expressed as stated in 3.22:

(7.6)

And taking into account that ${\textstyle \mathrm {W} _{z}}$ results from the contribution of the torque generated by the engine and the elastic reaction torque induced by the structure:

(7.7)

And recall how the elastic reaction torque was computed:

(7.8)

Accounting only for the torque in the rotation axis, and including equation 3.21 equation 7.2, we arrive at:

{\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} }})\\\\\mathrm {J} _{z}\;\alpha ={\frac {\textrm {P}}{\Omega }}-\sum _{j=1}^{n_{pt}}\left({\begin{bmatrix}\Delta {\widetilde {\mathbf {r} }}\\\vdots \end{bmatrix}}\left[\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} }\right]\right)_{\,(3\cdot j)}\end{matrix}}

(7.9)

The main obstacle in the resolution of the above system is that all ${\textstyle \mathbf {K} }$ , ${\textstyle \mathbf {D} }$ and ${\textstyle \mathbf {F} }$ depend on both ${\textstyle \Omega }$ and ${\textstyle \alpha }$ . Moreover, the way in which the previous matrices are defined (see 3.13) make the equations highly nonlinear in ${\textstyle \Omega }$ and ${\textstyle \alpha }$ . Hence, a proper numerical method has to be proposed in order to tackle the presented system of non-linear differential equations.

7.3 Numeric integration

In the context of coupled rigid-elastic dynamics, the unknowns are the vectors ${\textstyle \mathbf {d} }$ , ${\textstyle {\dot {\mathbf {d} }}}$ and ${\textstyle {\ddot {\mathbf {d} }}}$ regarding the vibration motion, and the scalars ${\textstyle {\theta }}$ , ${\textstyle {\dot {\theta }}}$ and ${\textstyle {\ddot {\theta }}}$ regarding the rotation motion. Using the presented Newmark ${\textstyle {\boldsymbol {\beta }}}$ -method, the previous system of differential equation can be transformed into a system of algebraic equations. The transformation this method proposes has already been reviewed in 5.2 for the elastic variables, and the only thing left is to apply it to the angular variables:

(7.10)

Hence, what is left is a system of non-linear algebraic equations with the unknowns ${\textstyle \alpha ^{n+1}}$ and ${\textstyle {\mathbf {d} }^{n+1}}$ . This problem could be tackled through a non-linear solving algorithm such as Newton-Raphson, widely used in computational mechanics. It is an iterative procedure based on the solving of a sequence of linear models [12]. The Newton method can in fact be applied to arbitrary systems of algebraic equations with ${\textstyle n}$ unknowns. However, its implementation is not an easy task and requires a great understanding of algebra. Moreover, as it is an iterative algorithm, it requires the different matrices to be computed more than once at each timestep, increasing considerably the computation time.

A simpler approach is based on the explicit central difference scheme, which computes accelerations at ${\textstyle {\textrm {t}}^{n}}$ and uses this value to update velocities and displacements at ${\textstyle {\textrm {t}}^{n+1}}$ . Although very easy to program, explicit methods require from very little timesteps to ensure conditional stability. Thus, computation time will arise with respect to the previous case of known kinematics. The general algorithm can be summarised as:

Obtain ${\textstyle {\mathbf {d} }^{n+1}}$ solving the system of equations stated in 5.34 (using ${\textstyle \Omega ^{n}}$ and ${\textstyle \alpha ^{n}}$ ) .
Once displacements at ${\textstyle {\textrm {t}}^{n+1}}$ are known, compute the elastic reaction torque ${\textstyle \mathrm {W} _{el}^{n+1}}$ from 7.9.
Find ${\textstyle \alpha ^{n+1}}$ from equation 7.7 and update ${\textstyle \Omega ^{n+1}}$ using expression 7.10.

The initial value of angular accelerations is found by solving the non-linear algebraic equation that involves initial conditions ${\textstyle \mathbf {d} _{0}}$ , ${\textstyle {\dot {\mathbf {d} }}_{0}}$ , ${\textstyle {\ddot {\mathbf {d} }}_{0}}$ and ${\textstyle \Omega _{0}}$ . As ${\textstyle \alpha _{0}}$ is the only unknown, the equation is solved by using a simple fixed-point iteration, using as initial guess ${\textstyle \alpha _{0}=0}$ . It has been found that four iterations suffice to obtain a converged solution.

7.4 Aerodynamic forces

Performing a simulation that accounts for fluid-structure interaction is far out of the scope of this section. Instead, the aerodynamic forces acting on the blade will be modelled using simplified aerodynamic relations. As the rotor is modelled in two dimensions, no forces acting perpendicular to the plane will be taken into consideration. Following the hypotheses presented in 7.1, the only force acting on the rotor plane is the parasite drag.

This aerodynamic resistance force depends on many parameters, such as the angle of attack, the surface roughness of the blade, the Mach number and the Reynolds number — a non-dimensional parameter that measures the strength of inertial forces of the fluid with respect to viscous ones. In order simplify things up, the effect of the latter variables will be encapsulated in a non-dimensional drag coefficient, ${\textstyle c_{D}}$ , which will be considered known. Otherwise, this analysis would be impracticable. The drag differential ( ${\textstyle \mathrm {d} {\textrm {D}}}$ ) acting on a surface differential ( ${\textstyle \mathrm {d} {\textrm {A}}}$ ) of the blade is given by:

(7.11)

Where ${\textstyle {\textrm {v}}}$ stands for the relative velocity of the surface differential with respect the fluid. As the blades are rotating in a plane, their solid rigid velocity is proportional to ${\textstyle \Omega \,\mathbf {r} _{i}}$ , where ${\textstyle \mathbf {r} _{i}}$ stands for the initial position. Taking into account the velocity induced by the elastic vibrations, and the change in radial position due to elastic displacements:

\mathbf {v} ^{2}=({\widetilde {\Omega }}\,(\mathbf {r} _{i}+\mathbf {u} )+{\dot {\mathbf {u} }})^{2}=\left[{\dot {\mathbf {u} }}^{2}+2\,({\widetilde {\Omega }}\,\,\mathbf {u} )\,{\dot {\mathbf {u} }}\,+2\,({\widetilde {\Omega }}\,\mathbf {r} _{i})\,\,{\dot {\mathbf {u} }}\,+({\widetilde {\Omega }}\,\mathbf {u} )^{2}+2\,({\widetilde {\Omega }}\,\mathbf {r} _{i})\,({\widetilde {\Omega }}\,\mathbf {u} )\,+({\widetilde {\Omega }}\,\mathbf {r} _{i})^{2}\right]

(7.12)

As the geometry is 2D, the drag force can be modelled in terms of the FEM as a body force, proportional to the area of each element. In three-dimensional problems this approach would not be correct, as aerodynamic forces only act on the surface and hence are to be modelled as boundary tractions. With this in mind, the term ${\textstyle \mathrm {d} {\textrm {D}}}$ is to be inserted as a body force ${\textstyle {\boldsymbol {f}}}$ in equation 3.13. As a result, the following expression is obtained:

(7.13)

The problem is that from ${\textstyle {\boldsymbol {f}}_{D}}$ will not only emerge components proportional to ${\textstyle {\dot {\mathbf {u} }}}$ and ${\textstyle \mathbf {u} }$ , but also non-linear terms such as ${\textstyle {\dot {\mathbf {u} }}^{2}}$ , ${\textstyle \mathbf {u} ^{2}}$ and ${\textstyle {\dot {\mathbf {u} }}\,\mathbf {u} }$ . The responsible for these nonlinearities to appear is the term ${\textstyle \mathbf {V} ^{2}}$ defined in 7.12. To make the problem more attainable, it will be supposed that ${\textstyle {\widetilde {\Omega }}\,\mathbf {r} _{i}\gg {\dot {\mathbf {u} }}}$ and ${\textstyle \mathbf {r} _{i}\gg \mathbf {u} }$ . Moreover, air density and drag coefficient will be considered constant along the domain. These hypotheses allow ${\textstyle {\boldsymbol {f}}_{D}}$ to be simplified as:

(7.14)

The above integral is decomposed into elemental matrices as reviewed in B, and computed using Gauss Quadrature. The major drawback of this simplified scenario is that drag forces no longer account for the effect of vibration, and thus air does not damp oscillating motions. In order to overcome this limitation, internal Rayleigh damping will be taken into consideration (recall equation 3.15). This cushioning will induce a new sources of stresses as expected from equation 2.15:

(7.15)

To close down this section, let's make an estimation of how much power is required to compensate the parasite drag when operating at nominal speed, say 400 rpm. Considering a constant chord ${\textstyle c}$ , the power the aerodynamic resistance of each blade induces is:

(7.16)

Considering flight at sea level ( ${\textstyle \rho _{\,{\hbox{air}}}=1.225\,Kg/m^{3}}$ ), a mean chord ${\textstyle c=0.45\,m}$ and a drag coefficient ${\textstyle c_{D}=0.008}$ , each blade requires from ${\textstyle 25.32\,kW}$ to surpass aerodynamic resistance. Nevertheless, the engine will need to inject more power to counteract the elastic reaction torque.

7.5 Rotor starting simulation

Once the FEM software is updated with the developed models from the previous sections, it is time to test it out simulating the starting of a rotor blade. Imagine that the studied helicopter has the following properties:

Number of blades: 3
Mean chord: 0.45
Nominal speed: 400 rpm
Motor nominal torque: 180 kN ${\textstyle \cdot }$ m
Engine nominal power: 580 kW
Nominal torque given at: ${\textstyle t_{a}}$ = 25 ms
Engine starting time: ${\textstyle t_{b}}$ = 50 ms
Nominal power given at: ${\textstyle t_{c}}$ = 120 ms

From basic helicopter theory, it is found that for hovering flight, a simple approximation of the power required to compensate the weight is:

(7.17)

Hence, considering an helicopter with a mass ${\textstyle m=2800\;Kg}$ , the power required to ensure hovering conditions at sea level is ${\textstyle {\textrm {P}}_{\,{\hbox{lift}}}=\;328.18\,kW}$ . Subtracting this value to the nominal power, it is found that the available power to counteract the elastic reaction torque is around ${\textstyle 84\;kW}$ per blade, which is higher than the parasite power induced by aerodynamic drag obtained in equation 7.16.

To be conservative, it will be supposed that the engine is not operating at nominal power but at its 60% ( ${\textstyle 50\;kW}$ per blade). Regarding the electric motor, it will be considered that at the beginning of the rotation motion the blades are undeformed, and hence all the torque — 60 kN ${\textstyle \cdot }$ m per blade — is used to accelerate the rotor. For sure, this development is an oversimplification of the real scenario, but still is good enough to provide a glimpse on how the blade dynamically behaves.

As the numerical integration requires from very small timesteps ( ${\textstyle \Delta {\textrm {t}}=5\mu \,s}$ ), the timespan is limited to ${\textstyle 5}$ s. To reduce computation time, stresses will only be computed every 20 timesteps. Dynamic displacements will also be saved at the same rate to reduce memory consumption. The integration is going to be performed using the constant average acceleration scheme. The considered structural Rayleigh damping is given by ${\textstyle {\overline {\alpha }}=0.25}$ and ${\textstyle {\overline {\beta }}=5\,\cdot \,10^{-6}}$ . With these inputs, the simulation takes up to three hours to complete.

The full results of the simulations are presented in the appendix (section D.3). The temporal evolution of the presented variables is closely linked with the way torque is injected into the blade, as angular accelerations are proportional to it. Different input torque will result in a completely different dynamic behaviour of the blade. For the studied case, torque and power are displayed in Figure 7.3.

Figure 7.3: Injected power and torque as a function of time.

Considering displacements at the tip of the blade as representative from the vibration motion, it is found that the blade is not undergoing continuous oscillation around and equilibrium configuration. Although it is true that oscillations are governing at the very beginning, displacements in ${\textstyle \varrho }$ tend to reach a non vibrating steady state. This effect was not sighted in previous simulations, where vibrations persisted in time. In all likelihood, the reduction of high frequency vibrations is induced by the structural damping acting on the blade. In fact, the SVD has been able to recover predominant damped oscillations — see RSVs 3 and 4 in Figure D.37 — that would otherwise be impossible to spot.

Much the same can be said about local axial displacements, but this time their tendency is not towards equilibrium but into a raise in time. This same behaviour has been spotted in earlier simulations, as the centrifugal force inducing axial displacements increases with ${\textstyle \Omega ^{2}}$ . Both displacements in ${\textstyle \varsigma }$ and angular speed are compared in Figure D.34.

Regarding the forces acting on the blade, there is a considerable peak at the beginning of the rotation motion, when the motor is injecting torque to the standstill rotor. During these first instants, all forces, accelerations and even stresses undergo oscillations. As time passes by and the engine achieves its nominal conditions, both reaction forces and accelerations stabilise into a considerably lower value. The only exception has to do with the torque induced by aerodynamic drag, that increases with the angular speed. In Figure D.35, angular acceleration is graphed along with the reaction torque, showing how acceleration increases when reactions decrease.

Figure 7.4: Deformed configuration (

) and

(KPa) of the blade for the moment of maximum stress.

When reviewing maximum stresses, it is found that the peaks correspond to the maximum deformed configuration. Again, this value reaches a steady condition as time passes by. What it has been found when carefully checking the results is that, after a few seconds, stresses tend to raise again. This effect is caused by the increase of aerodynamic drag, that induces forces proportional to ${\textstyle \Omega ^{2}}$ . This same condition is found in transversal displacements, hinting that there is a certain value of angular speed at which stresses and deformations are minimum. This effect can be spotted in Figure 7.5, where these minimums are represented with a red cross.

Figure 7.5: Temporal evolution of

displacements at the tip of the beam and maximum

.

It is worth pointing that in the simulated timespan, the angular velocity of the rotor has not reached its nominal value, but further timespans are far too restrictive in terms of computation time. The rotor would eventually surpass the nominal speed, as no control mechanism has been considered. To overcome this issue, injected power needs to be reduced when approaching the desired conditions.

8. Closure

Rotordynamics has been a field of great interest for centuries, being credited for the huge improvement in rotating machinery since the industrial revolution. Modern rotating components such as turbines, rotor blades and engines are still being designed using classical models developed centuries ago. What really created a game-changing shift in rotordynamics was the emergence of numerical methods in the earlier sixties. Implementation of the FEM gave rise to a completely new approach for structural analysis, allowing complex geometries to be solved in very short amounts of time.

Nowadays, the use of numerical software in structural analysis is a quite common solution. Despite of the great improvement in numerical capabilities over the last years, the FEM still requires from great processing power to perform complex simulations. Hence, it is not surprising that the actual tendency goes in the direction of more efficient methods instead of more accurate models. In this context, the present report has reviewed some of the available numerical tools to tackle the equations describing the vibration behaviour of a rotatory blade.

This thesis started with the elastic modelling of a simplified rotor blade undergoing one great rotation. The model supposes small strains and neglects the gyroscopic effect. When the obtained equations are posed in terms of the FEM, a softening effect is unfolded. The latter effect causes natural frequencies to decrease with rotation speed. However, real rotors do not exhibit this behaviour, but instead are affected by an opposite stiffening effect. Hence, the simplified model is not capable to capture the complete vibration phenomena, but even so it is a great tool that offers a glimpse on how rotordynamics solvers work.

For the static case, the developed FEM tool has proven to provide correct results, as it has been validated with analytic expressions. Regarding the dynamic scenario, the present work has focused on the numeric methods itself rather than on the description of the physical model. The goal of this thesis is not to bring to light new discoveries in the field of rotordynamics, but to give a glimpse on how numeric schemes work and in which way could they be more performing. The problem with structural analysis is that the governing equations rapidly become nonlinear, being numerical integration the only possible path to handle them.

For that purpose, a numerical scheme based on the Newmark- ${\textstyle \beta }$ method has been developed to study the dynamic response of rotating structures. It has been found that the rotation condition increases considerably the computational costs of the simulations, as matrices linked with rotating effects need to be computed and assembled at each timestep. In order to reduce simulation times, a more efficient assembly method has been used instead of the traditional approach, leading to an increase in performance by a factor of 15.

On the opposite side, when the physics describing the dynamic behaviour of a system are simple, modal decomposition analysis can be used to study its vibration behaviour. Modal analysis is a widespread tool that unfolds relevant information about the system with low computation costs. In structural analysis, modal decomposition is used to recover natural frequencies and mode shapes. The biggest drawback of this method is its inherent limitation to work only with simple dynamic systems.

Being aware of this situation, the truncated SVD has been found to be an alternative method that combines both numerical integration and modal decomposition to recover predominant patterns of the dynamic behaviour of the structure. If applied to structural analysis, it has been proven to recover predominant modes and frequencies of the system with great accuracy. The studied cases implied 2D and 3D rotating cantilever beams, although it could be extended to any arbitrary geometry. In order to recover oscillation frequencies from predominant patters, the developed method uses the discrete Fourier transform to catch the higher components in the frequency domain. One has, however, to keep in mind that the SVD is a probabilistic method, and thus gains reliability as the size of the sample increases.

The reader could come to think that the SVD entails no real advantage, as the input data matrix is still computed using numeric integration. Nothing further from the truth: the proposed method allows a deeper understanding of the vibration phenomena as it breaks up the whole problem into simpler components, easier to interpret and to postprocess. The postprocessing becomes thereby much simpler and faster, as the complete behaviour of the structure can be decomposed into a linear combination of predominant modes and frequencies. Hence, the problem numerical integration has about encountering multiple time and space scales in the postprocessing step is overcome. This make results analysis an easier and quicker task for the engineer, reducing both the time and computation resources consumed.

When testing the SVD in rotating beams, it has been found that the method is able to keep track of the softening effect predicted by the model. This capacity of capturing the reduction in vibration frequencies as angular velocity increases would be unthinkable through traditional approaches. Moreover, it has been spotted that the dynamic displacements can be approximated using only a few singular values, allowing an extensive reduction in postprocessing requirements.

In order to locate the developed work in the context of aerospace engineering, the built software has been adapted to work together with a rigid body model. This has allowed the transient study of the starting of an helicopter's rotor, giving plausible results and proving the actual capabilities of the method. When thinking about the potential capabilities of the SVD in structural analysis, it is too soon to reach a conclusion. The work is still in a very early state, and much investigation needs to be carried in this subject. This report has only accounted for a very simple rotor model, so nothing can be said about the validity of the results as there is no way they could be compared with empirical data. Until the proposed method is further developed, the present thesis will remain as a theoretical but plausible approach for real rotordynamics studies.

A. Theoretical frame

A.1 Introduction to elasticity

A.1.1 The elastic problem

Fundamentals of Theory of Elasticity date from the seventeenth century, back when Robert Hooke posed the law of elasticity, which states that the stretching of a solid is proportional to the force applied to it. This law was latter generalised to three dimensional bodies by Cauchy, introducing six components of stress ${\textstyle \sigma }$ , linearly related to six components of strain ${\textstyle \varepsilon }$ . Cauchy was the first to introduce the notion of stress at a point, given by the tractions per unit area across all plane elements trough the point.

This linear behaviour of deformable bodies is governed by 15 coupled partial differential equations, first derived by Navier. Even today, no closed form of the solution exists, even for simple structures. Elasticians have used energy principles to tackle this problem, which have become a fundamental tool to study elastic bodies, including those problems involving vibrations (see A.2). The Navier-Cauchy equations are [1]:

Equilibrium equations:

{\begin{aligned}{{\frac {\partial \sigma _{xx}}{\partial x}}+{\frac {\partial \sigma _{yx}}{\partial y}}+{\frac {\partial \sigma _{zx}}{\partial z}}+b_{x}=0}\\\\{{\frac {\partial \sigma _{xy}}{\partial x}}+{\frac {\partial \sigma _{yy}}{\partial y}}+{\frac {\partial \sigma _{zy}}{\partial z}}+b_{y}=0}\\\\{{\frac {\partial \sigma _{xz}}{\partial x}}+{\frac {\partial \sigma _{yz}}{\partial y}}+{\frac {\partial \sigma _{zz}}{\partial z}}+b_{z}=0}\end{aligned}}

(A.1)

Strain displacement relations:

{\begin{aligned}\varepsilon _{xx}&={\frac {\partial u}{\partial x}},\quad \varepsilon _{yz}={\frac {1}{2}}\left({\frac {\partial w}{\partial y}}+{\frac {\partial v}{\partial z}}\right)=\varepsilon _{zy}\\\\\varepsilon _{yy}&={\frac {\partial v}{\partial y}},\quad \varepsilon _{zx}={\frac {1}{2}}\left({\frac {\partial u}{\partial z}}+{\frac {\partial w}{\partial x}}\right)=\varepsilon _{xz}\\\\\varepsilon _{zz}&={\frac {\partial w}{\partial z}},\quad \varepsilon _{xy}={\frac {1}{2}}\left({\frac {\partial v}{\partial x}}+{\frac {\partial u}{\partial y}}\right)=\varepsilon _{yx}\end{aligned}}

(A.2)

Compatibility relations:

{\begin{aligned}{\frac {\partial ^{2}\varepsilon _{xx}}{\partial y^{2}}}+{\frac {\partial ^{2}\varepsilon _{yy}}{\partial x^{2}}}=2{\frac {\partial ^{2}\varepsilon _{xy}}{\partial x\partial y}},\quad {\frac {\partial ^{2}\varepsilon _{xx}}{\partial y\partial z}}={\frac {\partial }{\partial x}}\left(-{\frac {\partial \varepsilon _{yz}}{\partial x}}+{\frac {\partial \varepsilon _{zx}}{\partial y}}+{\frac {\partial \varepsilon _{xy}}{\partial z}}\right)\\\\{\frac {\partial ^{2}\varepsilon _{yy}}{\partial z^{2}}}+{\frac {\partial ^{2}\varepsilon _{zz}}{\partial y^{2}}}=2{\frac {\partial ^{2}\varepsilon _{yz}}{\partial y\partial z}},\quad {\frac {\partial ^{2}\varepsilon _{yy}}{\partial z\partial x}}={\frac {\partial }{\partial y}}\left({\frac {\partial \varepsilon _{yz}}{\partial x}}-{\frac {\partial \varepsilon _{zx}}{\partial y}}+{\frac {\partial \varepsilon _{xy}}{\partial z}}\right)\\\\{\frac {\partial ^{2}\varepsilon _{zz}}{\partial x^{2}}}+{\frac {\partial ^{2}\varepsilon _{xx}}{\partial z^{2}}}=2{\frac {\partial ^{2}\varepsilon _{zx}}{\partial z\partial x}},\quad {\frac {\partial ^{2}\varepsilon _{zz}}{\partial x\partial y}}={\frac {\partial }{\partial z}}\left({\frac {\partial \varepsilon _{yz}}{\partial x}}+{\frac {\partial \varepsilon _{zx}}{\partial y}}-{\frac {\partial \varepsilon _{xy}}{\partial z}}\right)\end{aligned}}

(A.3)

Equilibrium equations are obtained considering a general state of stress at an arbitrary point in the deformable body through an infinitesimal approach (Figure A.1). Note that shear stresses are represented with a ${\textstyle \tau }$ . If one approximates the spatial variations of the stresses in terms of first-order Taylor series, the resulting equilibrium equations are the ones posed in equation A.1. The strain-displacement equations link both strain and displacement field, defining the second as the derivative of the first. However, this is not enough, as another condition needs to be met: the body is still continuous after the deformation. The latter is ensured by the compatibility relations, which involve second derivatives.

Figure A.1: General stress state for a 3D deformable body [16]

A.1.2 Generalised Hooke's law

Robert Hooke was one of the first scientists to notice elastic behaviour of materials. He studied springs under traction and compression, and tried to pose an equation linking both the applied force and the elongation of the spring. He discovered that, for small deformations ${\textstyle x}$ , the required force was ${\textstyle F=kx}$ , where ${\textstyle k}$ is a constant characteristic of the spring. This law was later generalised to describe elastic objects in general, posing that the strain ${\textstyle \varepsilon }$ (that is, nondimensional deformation) is proportional to the stress ${\textstyle \sigma }$ (that is, force per unit area) applied to it.

Remembering the equations from the previous section, one can realise that general stresses and strains have multiple independent components, and thus the proportionally factor is no longer a single value, but a tensor. The tensorial form of Hooke's law introduces a concept that might not be intuitive: some elastic bodies deform in one direction when subjected to a force with a different direction. In such cases, proportionality will hold: if the directions are not changed, the deformation will increase proportional to the applied force. For that cases, the formulation in three dimensions has the following vectorial form:

\left[{\begin{array}{c}{F_{1}}\\{F_{2}}\\{F_{3}}\end{array}}\right]=\left[{\begin{array}{ccc}{k_{11}}&{k_{12}}&{k_{13}}\\{k_{21}}&{k_{22}}&{k_{23}}\\{k_{31}}&{k_{32}}&{k_{33}}\end{array}}\right]\left[{\begin{array}{c}{x_{1}}\\{x_{2}}\\{x_{3}}\end{array}}\right]

(A.4)

When moving from a discrete spring to continuous media, the formulation changes a little as the strain and stress state around a point cannot longer be described by a vector as all pushing, pulling and shearing effects have to be taken into account. This complexity is captured using a second order tensor, and in the tensorial formulation, forces and displacements are replaced by stresses and strains, respectively. In a Cartesian base, stresses and strains are represented with a first order tensor each:

{\boldsymbol {\varepsilon }}=\left[{\begin{array}{lll}{\varepsilon _{xx}}&{\varepsilon _{xy}}&{\varepsilon _{xz}}\\{\varepsilon _{yx}}&{\varepsilon _{yy}}&{\varepsilon _{yz}}\\{\varepsilon _{zx}}&{\varepsilon _{zy}}&{\varepsilon _{zz}}\end{array}}\right];\quad {\boldsymbol {\sigma }}=\left[{\begin{array}{lll}{\sigma _{xx}}&{\tau _{xy}}&{\tau _{xz}}\\{\tau _{yx}}&{\sigma _{yy}}&{\tau _{yz}}\\{\tau _{zx}}&{\tau _{zy}}&{\sigma _{zz}}\end{array}}\right]

(A.5)

Stress and strain tensors vary from a point to another inside continuum materials: ${\textstyle \mathbf {\varepsilon } }$ gives information about the displacements in the neighbourhood of the point, while ${\textstyle \mathbf {\sigma } }$ stands for the forces per unit area that neighbouring parcels exert on each other. The tensor linking ${\textstyle \mathbf {\varepsilon } }$ and ${\textstyle \mathbf {\sigma } }$ is known as the elasticity tensor ${\textstyle \mathbf {C} }$ , represented by a matrix of 81 real numbers ${\textstyle C_{ijkl}}$ :

(A.6)

Due to inherent symmetries of the model, the elasticity tensor ${\textstyle \mathbf {C} }$ can be reduced to 21 independent coefficients, while ${\textstyle \mathbf {\varepsilon } }$ and ${\textstyle \mathbf {\sigma } }$ are found to be symmetrical matrices ( ${\textstyle \sigma _{ij}=\sigma _{ji}}$ and ${\textstyle \varepsilon _{kl}=\varepsilon _{lk}}$ ), having thus only 6 independent values. In this context, the Voigt notation is introduced: strain and stress tensors are substituted by vectors including only independent coefficients. For the case of shear strains, they are substituted by engineering strains ${\textstyle \gamma }$ in Voigt notation ( ${\textstyle 2\varepsilon _{kl}=\gamma _{kl}}$ ) :

{\boldsymbol {\sigma }}\Rightarrow \quad \left[{\begin{array}{l}{\sigma _{x}}\\{\sigma _{y}}\\{\sigma _{z}}\\{\tau _{yz}}\\{\tau _{xz}}\\{\tau _{xy}}\end{array}}\right];\quad \quad {\boldsymbol {\varepsilon }}=\nabla ^{s}\mathbf {\underline {u}} \Rightarrow \left[{\begin{array}{c}{\varepsilon _{x}}\\{\varepsilon _{y}}\\{\varepsilon _{z}}\\{\gamma _{yz}}\\{\gamma _{xz}}\\{\gamma _{xy}}\end{array}}\right]=\left[{\begin{array}{c}{\frac {\partial u_{x}}{\partial x}}\\{\frac {\partial u_{y}}{\partial y}}\\{\frac {\partial u_{z}}{\partial z}}\\{{\frac {\partial u_{y}}{\partial z}}+{\frac {\partial u_{z}}{\partial y}}}\\{{\frac {\partial u_{x}}{\partial z}}+{\frac {\partial u_{x}}{\partial y}}}\\{{\frac {\partial u_{y}}{\partial x}}+{\frac {\partial u_{x}}{\partial y}}}\end{array}}\right]

(A.7)

In the above equation, ${\textstyle \nabla ^{s}\mathbf {\underline {u}} }$ stands for the symmetric gradient operator, which in 3D Cartesian coordinates adopts the following form:

\nabla ^{s}:=\left[{\begin{array}{lll}{\frac {\partial }{\partial x}}&{0}&{0}\\{0}&{\frac {\partial }{\partial y}}&{0}\\{0}&{0}&{\frac {\partial }{\partial z}}\\{0}&{\frac {\partial }{\partial z}}&{\frac {\partial }{\partial y}}\\{\frac {\partial }{\partial z}}&{0}&{\frac {\partial }{\partial x}}\\{\frac {\partial }{\partial y}}&{\frac {\partial }{\partial x}}&{0}\end{array}}\right]

(A.8)

For the case of isotropic materials —that is, uniform properties in all directions—, it can be found that the elasticity matrix depends only on two coefficients, the Young's Modulus ${\textstyle E}$ (standing for the material stiffness or resistance to be elastically deformed) and the Poisson's coefficient ${\textstyle \nu }$ . The latter is the negative of the ratio of transverse strain (unit lateral contraction) to axial strain, and the responsible for the Possion effect: the phenomenon in which a given material tends to expand in directions perpendicular to the direction of compression, and to contract in the directions transverse to the direction of stretching [16].

The Young's Modulus has units of pressure, while the Poisson's coefficient is dimensionless, and theoretically limited to the interval ${\textstyle -1<\nu {<0},5}$ , although in reality only positive values of ${\textstyle \nu }$ are found. However, ${\textstyle \mathbf {C} }$ is commonly expressed as a function of the Lame parameters ${\textstyle \lambda }$ and ${\textstyle \mu }$ . In indicial form:

(A.9)

In Voigt notation and for isotropic materials, Hooke's law adopts the form

\left[{\begin{array}{l}{\sigma _{x}}\\{\sigma _{y}}\\{\sigma _{z}}\\{\tau _{yz}}\\{\tau _{xz}}\\{\tau _{xy}}\end{array}}\right]=\left[{\begin{array}{cccccc}{2\mu {+\lambda }}&{\lambda }&{\lambda }&{0}&{0}&{0}\\{\lambda }&{2\mu {+\lambda }}&{\lambda }&{0}&{0}&{0}\\{\lambda }&{\lambda }&{2\mu {+\lambda }}&{0}&{0}&{0}\\{0}&{0}&{0}&{\mu }&{0}&{0}\\{0}&{0}&{0}&{0}&{\mu }&{0}\\{0}&{0}&{0}&{0}&{0}&{\mu }\end{array}}\right]\left[{\begin{array}{c}{\varepsilon _{x}}\\{\varepsilon _{y}}\\{\varepsilon _{z}}\\{\gamma _{yz}}\\{\gamma _{xz}}\\{\gamma _{xy}}\end{array}}\right]

(A.10)

The same can be deduced analogously for orthotropic and even anisotropic materials, with the difference that the constants ${\textstyle E}$ and ${\textstyle \nu }$ are no longer unique, but have a different value for each set of directions. Even thermal effects and other phenomena can be added into the equations describing the model. However, although Hooke's law is widely used, it presents great limitations as nonlinearities emerge when the deformations are big enough to be no longer considered small. Thus, nonlinear models have to be used, and the phenomena of plasticity has to be taken into account.

A.2 Classical analysis methods

The study of the vibration behaviour of rotatory blades is a well-known problem and the main subject of rotordynamics, a specialised branch of applied mechanics. In this section, a review on the main methods and models used over the last century for the analysis of rotating structures will be addressed, following the key concepts presented in [1].

A.2.1 Energy methods

One of the most fundamentals principles in Physics states that energy in the Universe is conserved. Energy can adopt many forms, but in the study of rotordynamics the most relevant are kinetic and potential (strain) energy. If energy changes from potential to kinetic repetitively in time, then vibration appears. In this section, the main methods used to describe the beam behaviour from an energetic point of view will be briefly presented.

Euler-Lagrange equations

The following equations take advantage of variational calculus to illustrate the behaviour of a one dimensional idealization of a cantilever beam. Let's consider a function that would represent the beam problem:

(A.11)

with ${\textstyle x}$ as independent variable and ${\textstyle y}$ and its derivatives as an admissible path. Considering the varied path defined by ${\textstyle {\breve {y}}=y+\varepsilon \eta }$ , being ${\textstyle \epsilon }$ a small parameter and ${\textstyle \eta }$ a differentiable function zero valued at ${\textstyle x_{1}}$ and ${\textstyle x_{2}}$ , equation A.11 over the varied path reads as:

(A.12)

If the function above is expanded:

(A.13)

Neglecting high order terms, one can state that ${\textstyle \left({\frac {\partial {\breve {I}}}{\partial \varepsilon }}\right)_{\varepsilon =0}=0}$ . Expanding and integrating by parts, the equation admits two solutions: or either ${\textstyle \eta }$ and ${\textstyle \eta ^{\prime }}$ are zero at ${\textstyle x_{1}}$ and ${\textstyle x_{2}}$ (Euler-Lagrange equation, eq. A.14) or its value is prescribed (boundary conditions, eq. A.15).

(A.14)

(A.15)

These equations describe one of the simplest structures, a one dimensional idealization of a cantilever beam. The greatness of this equation lies in converting the complex problem of solving 15 coupled differential equations (the Elasticity Problem) into a simple set of equations which are precise enough for engineering applications. For this particular case, the function ${\textstyle F}$ is found by computing the potential energy of the beam after a load is applied, and reads as:

(A.16)

Including eq. A.16 into eq. A.14 yields to

(A.17)

Better known as the beam equation: ${\textstyle w}$ states for lateral deflection, ${\textstyle x}$ for the independent coordinate along the beam axis, ${\textstyle E}$ for the Young's modulus, ${\textstyle I}$ for the second moment of area of the cross-section about the yy axis and ${\textstyle q}$ for a distributed load. An advantage of using energy methods instead of the Newtonian approach is that boundary conditions are completely fixed by equation A.15, whereas with classic methods they can be difficult to visualise.

Other methods

There are plenty of methods that use conservation of energy as basic principle to describe the behaviour of beams. The Lagrange method, for example, assumes that the solution of the equation satisfies the kinematic boundary conditions. For cantilever boundary conditions, in fact, a fourth grade polynomial is found to satisfy boundary conditions in (A.15). The solution is then found by minimizing total potential energy.

Lagrange also studied time-dependent structures, that is, the phenomena of vibration. The conservation of energy for a dynamic system states as ${\textstyle {\frac {d}{dt}}\left({\frac {\partial T}{\partial {\dot {q}}}}\right)-{\frac {\partial T}{\partial q}}+{\frac {\partial U}{\partial q}}=0}$ , ${\textstyle T}$ representing the kinetic energy and ${\textstyle U}$ the strain energy. If one uses ${\textstyle n}$ terms to define both kinetic and strain energy, the following matrix system is obtained, which leads to an eigenvalue problem that was not quite practical to solve until the ${\textstyle 20^{th}}$ century:

(A.18)

Rayleigh's approach, on the other hand, consisted in finding a method for calculating fundamental natural frequencies with a formulation that could be tabulated and used in industry as a routine numerical operations, so the designs could be faster developed. To do so, Rayleigh used the conservation principle which states that, for a vibrating conservative system in simple harmonic motion, the maximum potential energy at maximum displacement is equal to the maximum kinetic energy when the system passes trough its equilibrium with maximum velocity.

Another great contributor to this field was Grigoryevich Galerkin, who in 1915 published an article on how to solve differential equations boundary problems with an approximate solution method. Differential equations of the type of equation A.14 cannot be solved except in a very few special cases. The key here is to multiply equation A.14 by ${\textstyle \eta (x)}$ , which is assumed to be the solution itself, and then integrating this product from ${\textstyle x_{1}}$ to ${\textstyle x_{2}}$ . The solution is found by making this integral zero. In the selection of ${\textstyle \eta (x)}$ , the Galerkin method uses ${\textstyle {\breve {y}}}$ itself. As in the Lagrange method, if one aims to solve the beam equation, it will be necessary to approximate the deflection first.

Modern Finite Element Methods rely, however, on Hamilton's Principle, the most general form for all vibrating systems, and one of the most powerful when applied to rotordynamics. A general equation for beam rotordynamics is found in terms of this principle for constant angular acceleration ${\textstyle \alpha }$ :

{\begin{array}{c}{EI_{xx}y^{\prime \prime \prime \prime }+\rho A\left[{\begin{array}{c}{{\ddot {y}}-\left(\omega _{0}+\alpha t\right)^{2}\left\{y+R_{2}y^{\prime \prime }-(R+z)y^{\prime }\right\}+}\\{+2\left(\omega _{0}+\alpha t\right)\left\{{\dot {y}}y^{\prime }-y^{\prime \prime }\int _{z}^{l}{\dot {y}}dz-\int _{0}^{z}y^{\prime }{\dot {y}}^{\prime }dz\right\}}\\{+\alpha \left\{yy^{\prime }-{\frac {1}{2}}\int _{0}^{z}y^{\prime 2}dz-y^{\prime \prime }\int _{z}^{l}ydz\right\}}\end{array}}\right]-}\\\\{-\left(\rho I_{xx}+{\frac {k\rho EI_{xx}}{G}}\right){\ddot {y}}^{\prime \prime }+I_{xx}k\rho ^{2}{\ddot {y}}=-\rho A\alpha (R+z)}\end{array}}

(A.19)

Neglecting the effect of shear ( ${\textstyle k=0}$ for long blades, as in helicopters), the equation can be simplified and some terms identified (Figure A.2). The following equation is solved using approximate methods, and the Galerkin method is a powerful tool to put it in terms of Finite Element Method. As can be seen, Coriolis forces and angular acceleration add non-linearities to the equation. A stiffness stiffening and softening due to rotation can also be spotted.

Figure A.2: Simplified equation for long rotating beams [1]

A.2.2 Graphical methods

At the beginning of the ${\textstyle 20^{th}}$ , a rapid growth of rotating machinery made the need of rapid designs in the industrial environment arise. To handle this problem, graphical methods were used to determine natural frequencies of rotating structures. The goal was to create methods to be used by semi-skilled engineers and handled in the shortest possible time.

At that time, beam theory was already developed and used in many fields. A very common method to handle vibration problems consisted on drawing both the shear force and bending moment diagrams. With those, the deflection diagram could also be drawn and the integration was converted into graphical summation. This is a graphical way to interpret the Stodola-Viannello method, which considered deflection itself as a proportional load that produces the deflection. However, this process of drawing the diagrams is rather slow and depends on the accuracy of the drawings. That is why the Stodola-Viannello method became popular in the industrial world not in its graphical form, but as a tabulated iterative process.

In 1922, Holzer proposed another tabulated iterative method, based on Stodola's, which was able to capture higher torsional modes and its frequency. This method was widely used, and for instance a Railway Diesel-Electric Drive was designed. The main drawback of the Holzer method was the impossibility to predict deflection, slope, shear and bending moments.

The first graphical method used in aircraft wings and turbine blades was Myklestad’s method, widely spread during World War II. This tabular iterative method follows Holzer's in the way of capturing high torsional modes, but also computes deflection, slope, shear and bending moments. Phrol's method is quite similar to Myklestad’s, but addressed to shafts instead of aircraft structures. The main advantage of the latter was that the analysis included the rotatory inertia of the disks, thus adding the so-called gyroscopic effect.

This methods were widely used until the first half of the past century. With the advanced computational facilities, finite element methods gradually overtook these tabular and energy methods.

A.2.3 Matrix methods

As it has already been addressed in equation A.18, the Lagrange method enables the transformation of the differential equations governing the dynamics of free vibration systems into homogeneous algebraic equations. This leads to a well known eigenvalue problem, were the key is to find the vibration mode shapes (eigenvectors) and the natural frequencies (eigenvalues). However, the computational needs for practical systems were not present until the late ${\textstyle 20^{th}}$ century.

However, some methods were developed to handle these kind of systems with relative ease: Jacobsen and Ayre's, Gräffe's, Prieb's and Hahn's are just a few examples of methods that dealt with matrix systems to solve structures.

As matrix methods proven to be a good promise for vibration problems, Duncan and Collar developed an iterative method (essentially an extension of Stodola-Vianello's) that enabled to handle relatively large matrices. Transfer matrix methods also gained popularity because they were easier to adapt for computerization, as they only involve simple matrix multiplications, and the overall transfer matrix is always small. The latter method was extensively used in rotor dynamics calculations for determining critical speeds and unbalanced response. This method was extended even for transient whirl analysis.

But as processing power of computers started to grow, tabular and matrix methods lost their place and nowadays have no use in modern designs.

B. Elemental matrix computation

When the physical system of study is posed in terms of the FEM, spatial coordinates are discretised and the resulting equations are ruled by a set of global matrices (see in fact equation 3.14). That is why matrix computation is always an important step of every FEM solver. However, those matrices involve integrals over the domain which have no closed solution unless they are split into local or elemental coordinates. Under those conditions, the global matrices are decomposed into smaller, elemental matrices. If mapped into the parent domain, the integrals can be approximated using the Gauss Quadrature method. The aim of this chapter is to review the process of assembly: how global matrices are built from elemental ones. The approximation of the latter will also be assessed.

B.1 The static stiffness matrix

From equation 3.13, the expression of the stiffness matrix for static structures ${\textstyle \mathbf {K} _{\hbox{ static}}}$ is, in terms of the FEM:

(B.1)

The above integral is to be evaluated over the domain ${\textstyle \mho }$ , which becomes infeasible unless local formulation is adopted. Using the formulation presented in section 3.26, the above integral can be splitted into elemental integrals as follows:

(B.2)

Where ${\textstyle \mathbf {L} ^{e}}$ is a ${\textstyle n_{d}n^{e}\times n_{d}n_{pt}}$ matrix known as boolean connectivity matrix, consisting of the integers 0 and 1 relating element nodal variables ${\textstyle {\boldsymbol {d}}^{e}\in \mathbb {R} ^{n_{d}n^{e}}}$ with the global nodal vector ( ${\textstyle \mathbf {d} ^{e}=\mathbf {L} ^{e}\mathbf {d} }$ ). With all, the problem of computing a single ${\textstyle {n}_{d}{n}_{pt}\times {n}_{d}{n}_{pt}}$ global matrix reduces to the computation of ${\textstyle \mathbf {n} _{\mathbf {el} }}$ ${\textstyle {n}_{d}{n}^{e}\times {n}_{d}{n}^{e}}$ elemental matrices. For the sake of simplicity, the assembly operation is expressed as:

(B.3)

The only thing left is finding a reliable method to compute ${\textstyle \mathbf {K} _{\hbox{ static}}^{e}}$ . As integration boundaries may differ from one element to another, the best method to evaluate elemental matrices is by transforming them into the parent domain, where their geometry is normalised and thus Gauss Quadrature can be applied in the same way for all elements. Using the formulation presented in 3.3.1:

\mathbf {K} _{\hbox{ static}}^{e}=\int _{\mho }^{}\mathbf {B} ^{e^{T}}\mathbf {C} \;\mathbf {B} ^{e}\;d\mho ^{e}=\int _{\mho }^{}J^{e}\;\mathbf {B} ^{e^{T}}\mathbf {C} \;\mathbf {B} ^{e}\;d\mho _{\xi }=\sum _{g=1}^{m}{w}_{g}\left(J^{e}\mathbf {B} ^{e^{T}}\mathbf {CB} ^{e}\right)_{{\boldsymbol {\xi }}={\boldsymbol {\xi }}_{g}}

(B.4)

The number of Gauss points ${\textstyle m}$ required for the integration depends on the order of the employed shape functions, and thus is to be specified for each type of element.

B.2 The force vector

In equation 4.1, ${\textstyle \mathbf {F} }$ stands for the force vector accounting for the effect of the body forces ${\textstyle \mathbf {F} _{b}}$ , boundary tractions ${\textstyle \mathbf {F} _{t}}$ and centrifugal forces ${\textstyle \mathbf {F} _{\hbox{ rot}}}$ . The methodology used to compute them is the same as in the previous case: splitting the global vector into elemental vectors and mapping them into the parent domain, where Gauss Quadrature can be applied.

B.2.1 Body forces

Recall from equation 3.13 that body and rotation force vectors can be expressed in terms of the FEM as:

(B.5)

Where ${\textstyle {\boldsymbol {f}}}$ stands for the vector of external body forces acting on the structure. Notice that the shape of both ${\textstyle \mathbf {F} _{b}}$ and ${\textstyle \mathbf {F} _{\hbox{rot}}}$ is quite similar, and so they can be written as:

(B.6)

Being ${\textstyle {\widehat {\mathrm {f} }}}$ the total contribution of body and rotation forces. As in the previous case, the global vector can be splitted into elemental force vectors using the assembly operator:

(B.7)

The problem lies now in how to define the vector ${\textstyle {\widehat {\mathrm {f} }}}$ along each element. A common approach is to replace it by an approximation constructed by interpolation of its nodal values ${\textstyle {\widehat {\mathrm {f} }}\;^{e}}$ using the same shape functions as for the displacement field:

\left.{\widehat {\mathrm {f} }}\;\right|_{{\boldsymbol {x}}\in {\overline {\mho }}^{\;e}}\approx N^{e}\;\;{\widehat {\mathrm {f} }}\;^{e}\;,\;\;\;{\hbox{   with  }}\;\;\;{\widehat {\mathrm {f} }}\;^{e}=\left[{\begin{array}{c}{{\widehat {\mathrm {f} }}\;\left({\boldsymbol {x}}_{1}^{e}\right)}\\{{\widehat {\mathrm {f} }}\;\left({\boldsymbol {x}}_{2}^{e}\right)}\\{\vdots }\\{{\widehat {\mathrm {f} }}\;\left({\boldsymbol {x}}_{n^{e}}^{e}\right)}\end{array}}\right]

(B.8)

Note that this approximation can induce errors if, for example, two neighbouring elements have different densities, as density would be a nodal property rather than an elemental one. Using Gauss Quadrature, the sum of the elemental body and rotation force vectors is expressed as:

(B.9)

B.2.2 Traction forces

When determining the external forces due to boundary traction conditions ${\textstyle \mathbf {F} _{t}}$ , it is convenient to distinguish between point loads ${\textstyle \mathbf {F} _{\hbox{pnt}}}$ and distributed loads ${\textstyle \mathbf {F} _{\hbox{dis}}}$ . The determination of ${\textstyle \mathbf {F} _{\hbox{pnt}}}$ is straightforward: if a point force ${\textstyle \mathbf {Q} \in \mathbb {R} ^{n_{d}}}$ is applied to node ${\textstyle \mathbf {A} }$ , then ${\textstyle (\mathbf {F} _{\hbox{pnt}})_{A}=\mathbf {Q} }$ . Point loads are a particular case of distributed loads, where the force function is the Dirac delta operator: zero valued everywhere but A, where it presents a spike.

The study of distributed loads requires from a more meticulous approach, and a new formulation regarding boundary elements is to be posed. To start with, let's define ${\textstyle \mathbf {E} _{b}^{i}}$ as the set of elements lying on the Neumann boundary ${\textstyle \Gamma _{\sigma }^{i}}$ , with ${\textstyle i=1,\;\cdots \;,n_{d}}$ :

(B.10)

Being ${\textstyle n_{b}^{e}}$ the number of nodes per boundary element, and assuming for simplicity it is the same for all elements, one can define ${\textstyle \mathbf {CNb} ^{i}}$ as the boundary connectivity matrix, a ${\textstyle {n}_{el}^{b,i}\times {n}_{b}^{e}}$ matrix containing, in each row, the global indices of the boundary nodes of each patch ${\textstyle \Gamma _{\sigma }^{e,i}\left(e\in E_{b}^{i}\right)}$ . For each boundary element ${\textstyle \Gamma _{\sigma }^{e,i}}$ , a matrix of boundary shape functions is defined:

(B.11)

Notice that the boundary shape functions are different from the ones used to interpolate the displacements, as the dimension of the boundary is always smaller than the dimension of the elements, ${\textstyle {\overline {N}}_{a}^{e}:\Gamma _{\sigma }^{e,i}\rightarrow \mathbb {R} }$ . Recall the definition of the global vector of prescribed boundary tractions (equation 4.1 ), whose decomposition into elemental integrals reads as:

\mathbf {F} _{\hbox{dis}}=\sum _{i=1}^{n_{d}}\int _{\Gamma _{\sigma }^{i}}{\overline {\mathbf {N} }}_{i}^{T}\;{\overline {t}}^{\,i}d\Gamma =\sum _{i=1}^{n_{d}}\;\;\sum _{e\,\in \,\mathbf {E} _{b}^{i}}\left({\overline {\mathbf {L} }}^{e,i}\right)^{T}\int _{\Gamma _{\sigma }^{e,i}}{\overline {\mathbf {N} }}_{i}^{T}\;{\overline {t}}^{\,i}d\Gamma

(B.12)

Where ${\textstyle {\overline {\mathbf {L} }}^{e,i}}$ is the Boolean operator linking the global ${\textstyle \left({\mathbf {c} }\in \mathbb {R} ^{n_{d}\;n_{b}^{e}}\right)}$ and local ${\textstyle \left({\overline {\mathbf {c} }}^{e,i}\in \mathbb {R} ^{n_{b}^{e}}\right)}$ vectors of test coefficients. The above expression can be expressed as:

\mathbf {F} _{\hbox{dis}}=\sum _{i=1}^{n_{d}}\mathbf {F} _{\hbox{dis}}^{i}=\sum _{i=1}^{n_{d}}\;\;\sum _{e\,\in \,\mathbf {E} _{b}^{i}}\left({\overline {\mathbf {L} }}^{e,i}\right)^{T}\mathbf {F} _{dis}^{e,i}=\sum _{i=1}^{n_{d}}\;\;\sum _{e\,\in \,\mathbf {E} _{b}^{i}}\left({\overline {\mathbf {L} }}^{e,i}\right)^{T}\int _{\Gamma _{\sigma }^{e,i}}{\overline {\mathbf {N} }}_{i}^{T}\;{\overline {\mathbf {t} }}^{\,i}d\Gamma

(B.13)

Where ${\textstyle \mathbf {F} _{dis}^{e,i}\in \mathbb {R} ^{n_{b}^{e}}}$ is the elemental traction vector due to distributed loads along direction i, ${\textstyle \mathbf {F} _{\hbox{dis}}^{i}\in \mathbb {R} ^{n_{d}n_{pt}}}$ is the global vector due to distributed loads along direction i and ${\textstyle \mathbf {F} _{dis}\in \mathbb {R} ^{n_{d}n_{pt}}}$ is the global vector due to distributed loads. As for the case of body forces, the input function ${\textstyle {\overline {\mathbf {t} }}^{\,i}}$ is commonly replaced by an approximation constructed by interpolation of its nodal values ${\textstyle {\overline {\mathbf {t} }}^{\,e,i}}$ using the shape functions employed for interpolating the boundary test functions:

\left.{\overline {\mathrm {t} }}^{i}\;\right|_{{\boldsymbol {x}}\in \Gamma _{\sigma }^{\;e,i}}\approx {\overline {\mathbf {N} }}^{e}\;\;{\overline {\mathrm {t} }}^{\;e,i}\;,\;\;\;{\hbox{   with  }}\;\;\;{\overline {\mathrm {t} }}^{\;e,i}=\left[{\begin{array}{c}{{\overline {\mathrm {t} }}^{i}\;\left({\boldsymbol {x}}_{1}^{e}\right)}\\{{\overline {\mathrm {t} }}^{i}\;\left({\boldsymbol {x}}_{2}^{e}\right)}\\{\vdots }\\{{\overline {\mathrm {t} }}^{i}\;\left({\boldsymbol {x}}_{n_{b}^{e}}^{e}\right)}\end{array}}\right]

(B.14)

Using this interpolation, and taking advantage of the assembly operator, the global vector due to distributed loads along direction i results in:

(B.15)

If mapped into the parent domain, the computation of the elemental traction vector due to distributed loads along direction i can be approximated using the Gauss Quadrature method:

\mathbf {F} _{\hbox{dis}}^{e,i}=\int _{\Gamma _{\sigma }^{e,i}}{\overline {\mathbf {N} }}^{e^{T}}\left({\overline {\mathbf {N} }}^{e}{\overline {\mathbf {t} }}^{e,i}\right)d\Gamma =\int _{{\Gamma _{\sigma }}_{\xi }}J_{b}^{e,i}\;{\overline {\mathbf {N} }}^{e^{T}}\left({\overline {\mathbf {N} }}^{e}{\overline {\mathbf {t} }}^{e,i}\right)d\Gamma _{\xi }=\sum _{g=1}^{m_{b}}{\overline {w}}_{g}\left(J_{b}^{e,i}\;{\overline {\boldsymbol {N}}}^{e^{T}}\left({\overline {\boldsymbol {N}}}^{e}{\overline {\boldsymbol {t}}}^{e,i}\right)\right)_{{\overline {\xi }}={\overline {\xi }}g}

(B.16)

As it has already been stressed out, ${\textstyle \Gamma _{\sigma }^{e,i}\subset \mathbb {R} ^{n_{d}-1}}$ and thus the shape functions and weights appearing in the above expressions are different from those computed when assessing the stiffness matrix and the body forces vector. For 3D problems, the boundary elements are 2D, and therefore a different Gauss rule is applied.

B.3 The mass matrix

It is found from equation 3.13 that the mass matrix can be expressed in terms of the FEM as:

(B.17)

The process of splitting the above integral into elemental ones, and approximating the latter using the Gauss Quadrature method is identical as for the static stiffness matrix case. Using the assembly operator:

(B.18)

The elemental matrices are then mapped into the parent domain and approximated using Gauss Quadrature:

\mathbf {M} ^{e}=\int _{\mho }^{}\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\;d\mho ^{e}=\int _{\mho }^{}J^{e}\;\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\;d\mho _{\xi }=\sum _{g=1}^{m}{w}_{g}\left(J^{e}\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\right)_{{\boldsymbol {\xi }}={\boldsymbol {\xi }}_{g}}

(B.19)

B.4 The static damping matrix

It is found from equation 3.13 that the static damping matrix can be expressed in terms of the FEM as:

(B.20)

Using the approximation presented in equation 3.15, the above integral is transformed into:

(B.21)

Where ${\textstyle \mathbf {M} }$ and ${\textstyle \mathbf {K} \;_{\hbox{static}}}$ are known matrices, while ${\textstyle {\overline {\alpha }}}$ and ${\textstyle {\overline {\beta }}}$ are constant parameters.

B.5 The rotation matrices

Recall that the rotation of the structure induces new matrices to appear in equation 3.13:

(B.22)

Its computation method is not so different from the mass matrix case:

(B.23)

(B.24)

Using the Gauss Quadrature method to approximate the elemental integrals:

\mathbf {K} \;_{\hbox{rot}}^{e}=\int _{\mho }^{}J^{e}\;\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\mathbf {N} ^{e}\;d\mho _{\xi }=\sum _{g=1}^{m}{w}_{g}\left(J^{e}\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;({\widetilde {\Omega }}^{2}+{\widetilde {\alpha }})\;\mathbf {N} ^{e}\right)_{{\boldsymbol {\xi }}={\boldsymbol {\xi }}_{g}}

(B.25)

(B.26)

C. Example of vectorised assembly

The FEM is able to solve boundary value problems along very complex geometries because it splits the domain into a finite number of smaller partitions called elements. In numerical terms, this implies that global matrices are not computed all at once, but element by element. Then, all the elemental matrices need to be assembled, but traditional assembly algorithms can be very slow. For this reason a vectorised assembly algorithm has been used in the present project. The code was provided by prof. J. Hernández and adapted to the case of rotating structures. The goal of the present chapter is to illustrate how vectorised assembly works by developing, step by step, the assembly of the mass matrix for two linear quadrilateral elements.

C.1 Mass matrix vectorised assembly

To start with, recall the definition of the mass matrix and its decomposition into elemental components from the previous chapter:

{\begin{matrix}\mathbf {M} =\int _{\mho }^{}\mathbf {N} ^{T}{\boldsymbol {\rho }}\;\mathbf {N} \;d\mho \\\\\mathbf {M} ={\overset {n_{el}}{\underset {e=1}{\mathbf {A} }}}\mathbf {M} ^{e}\;\;,\;\;\;\;\mathbf {M} ^{e}=\int _{\mho }^{}\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\;d\mho ^{e}\\\\\mathbf {M} ^{e}=\int _{\mho }^{}\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\;d\mho ^{e}=\int _{\mho }^{}J^{e}\;\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\;d\mho _{\xi }=\sum _{g=1}^{m}{w}_{g}\left(J^{e}\mathbf {N} ^{e^{T}}{\boldsymbol {\rho }}\;\mathbf {N} ^{e}\right)_{{\boldsymbol {\xi }}={\boldsymbol {\xi }}_{g}}\end{matrix}}

(C.1)

To keep the example as simple as possible, let's consider the following geometry made from two linear quadrilateral elements (Figure C.1).

Figure C.1: 2D mesh composed by two quadrilateral elements.

Remember that for 2D elements, the numeration order is not important as well as it follows counter-clockwise direction. In the above figure, numbers inside a circle represent local numeration. When the geometric domain is not as trivial as the presented one, keeping track of the global numeration is not an easy task. On the other hand, local numeration is always performed in the same way. Hence, one may suspect that some properties remain undisturbed in the local frame.

Following the previous idea, the vectorised assembly takes advantage of the fact that in the parent domain, all the elements are geometrically identical and so are their shape functions. Hence, it is only necessary to compute the elemental matrix of one element and use it to construct a matrix including all the elemental shape functions matrices. The latter matrix, referred as ${\textstyle \mathbf {N} _{elem}}$ , is actually a ${\textstyle (n_{g}\;n_{d}\;n_{el}\times n^{e}\;n_{d})}$ matrix.

Remember that, for vector fields, the elemental matrix of shape functions ${\textstyle \mathbf {N} ^{e}}$ was defined as:

(C.2)

Where ${\textstyle \mathbf {I} }$ is the ${\textstyle n_{d}\times n_{d}}$ identity matrix and ${\textstyle N_{A}}$ is the scalar-valued shape function of node ${\textstyle A}$ . Let's define now ${\textstyle \mathbf {N} _{E}}$ , a ${\textstyle (n_{g}\;n_{d}\times n^{e}\;n_{d})}$ matrix that incorporates the value of ${\textstyle \mathbf {N} ^{e}}$ at every Gauss Point:

\mathbf {N} _{E}={\begin{bmatrix}\mathbf {N} _{\xi _{1}}^{1}&\mathbf {N} _{\xi _{1}}^{2}&\cdots &\mathbf {N} _{\xi _{1}}^{\,n_{e}}\\\\\mathbf {N} _{\xi _{2}}^{1}&\mathbf {N} _{\xi _{2}}^{2}&\cdots &\mathbf {N} _{\xi _{2}}^{\,n_{e}}\\\vdots &\vdots &\ddots &\vdots \\\mathbf {N} _{\xi _{n_{g}}}^{1}&\mathbf {N} _{\xi _{n_{g}}}^{2}&\cdots &\mathbf {N} _{\xi _{n_{g}}}^{\,n_{e}}\end{bmatrix}}={\begin{bmatrix}{\mathcal {}}{N}^{1}\\{\mathcal {}}{N}^{2}\\\vdots \\{\mathcal {}}{N}^{\,n_{pt}}\end{bmatrix}}

(C.3)

Where ${\textstyle {\mathcal {}}{N}}$ is a ${\textstyle n_{g}\;n_{d}\times n_{d}}$ matrix. Taking advantage of the fact that ${\textstyle \mathbf {N} _{E}}$ is the same for all elements, ${\textstyle \mathbf {N} _{elem}}$ can be easily constructed by repeating ${\textstyle \mathbf {N} _{E}}$ a number of ${\textstyle n_{el}}$ times:

(C.4)

The next step is to built ${\textstyle \mathbf {N} _{st}}$ , a ${\textstyle (n_{g}\;n_{d}\;n_{el}\times n_{pt}\;n_{d})}$ matrix that takes into account the global formulation and is constructed from ${\textstyle \mathbf {N} _{elem}}$ by using the connectivity matrix. For instance, for the presented geometry and with ${\textstyle n^{e}=4}$ and ${\textstyle n_{g}=4}$ , ${\textstyle \mathbf {N} _{st}}$ adopts the following form:

{\begin{matrix}\mathbf {N} _{st}={\begin{bmatrix}{\mathcal {}}{N}^{1}&{\mathcal {}}{N}^{2}&0&0&{\mathcal {}}{N}^{3}&{\mathcal {}}{N}^{4}\\0&{\mathcal {}}{N}^{1}&{\mathcal {}}{N}^{2}&{\mathcal {}}{N}^{3}&{\mathcal {}}{N}^{4}&0\end{bmatrix}}{\begin{matrix}e_{1}\\e_{2}\end{matrix}}\\{\begin{matrix}\;\;\;\;\;\;\;\;\;{\textrm {n}}_{1}&\;{\textrm {n}}_{2}&\;\;\;{\textrm {n}}_{3}&\,\;\;\;{\textrm {n}}_{4}&\;\;\;{\textrm {n}}_{5}&\,\;\;\;{\textrm {n}}_{6}\end{matrix}}\end{matrix}}

(C.5)

The grace of this approach is that once ${\textstyle \mathbf {N} _{E}}$ is computed, the transformations that follow are performed using optimised in-built functions such as sparse, reshape and repmat:

m = ndim*nelem*ngaus ; % Number of rows
n = nnode*ndim ;          % Number of columns
nzmaxLOC = m*nnodeE*ndim ;   % Maximum number of zeros
Nst = sparse([],[],[],m,n,nzmaxLOC); % Allocating memory for Bst

for inode = 1:nnodeE   % Loop over number of nodes at each element
    
    setnodes = CONNECT(:,inode) ;  % Global numbering of the
    %inode-th node of each element
      for idime = 1:ndim   % loop over spatial dimensions
        DOFloc = (inode-1)*ndim+idime ;  % Corresponding indices
        s = Nelems(:,DOFloc) ;% Corresponding column
        i = [1:m]' ;  % All rows
        DOFglo = (setnodes-1)*ndim+idime ; %  b) Columns
        j = repmat(DOFglo', ndim*ngaus,1) ;
        j = reshape(j,length(i),1) ;
        % Assembly process  (as sum of matrix with sparse representation)
        Nst = Nst + sparse(i,j,s,m,n,m) ;
    end
end

Until now, the Gauss weights and Jacobians were not considered. So let's define WSTs, a ${\textstyle n_{el}\times n_{g}}$ vector containing the product of weights and Jacobians at all Gauss points:

\mathbf {WSTs} ={\begin{bmatrix}s^{1}\\s^{2}\\\vdots \\{\hbox{s}}^{\,n_{el}}\end{bmatrix}}\;\;,\;\;s^{\,e}={\begin{bmatrix}w_{1}\;\cdot \;J_{\xi _{1}}^{e}\\w_{2}\;\cdot \;J_{\xi _{2}}^{e}\\\vdots \\w_{\,n_{g}}\;\cdot \;J_{\xi _{\,n_{g}}}^{e}\end{bmatrix}}

(C.6)

The previous matrix is to be transformed into diagonal. To do so, let's first define the following auxiliary matrices:

\mathbf {w} _{i}^{e}=\left(w_{i}\;\cdot \;\mathbf {J} _{\xi _{i}}^{e}\right)\mathbf {I} \;\;,\;\;\mathbf {W} _{d}^{e}={\begin{bmatrix}\mathbf {w} _{1}^{e}&[\,0\,]&[\,0\,]&[\,0\,]\\[][\,0\,]&\mathbf {w} _{2}^{e}&[\,0\,]&[\,0\,]\\\vdots &\vdots &\ddots &\vdots \\[][\,0\,]&\cdots &\cdots &\mathbf {w} _{n_{g}}^{e}&\end{bmatrix}}

(C.7)

And now the diagonalised version of WSTs reads as:

\mathbf {W} _{diag}={\begin{pmatrix}\mathbf {W} _{d}^{\,1}&[\,0\,]&\cdots &[\,0\,]\\\\[][\,0\,]&\mathbf {W} _{d}^{\,2}&\cdots &[\,0\,]\\\\\vdots &\vdots &\ddots &\vdots \\\\[][\,0\,]&[\,0\,]&\cdots &\mathbf {W} _{d}^{\,n_{el}}\end{pmatrix}}

(C.8)

To close down the section, let's take density into consideration. Recall that density is given as a vector input, having a value for each element:

(C.9)

The next step is to transform the previous expression into a global matrix, accounting for each dimension and Gauss point. To do so, let's consider the following matrices:

{\hbox{dens}}^{\,e}=\rho ^{\,e}\;\mathbf {I} \;\;,\;\;{{\hbox{dens}}_{\,\mathbf {G} }}^{\,e}={\overset {(n_{d}\;n_{d}\;n_{g})}{\overbrace {\begin{bmatrix}{\hbox{dens}}^{\,e}&[\,0\,]&\cdots &[\,0\,]\\[][\,0\,]&{\hbox{ dens}}^{\,e}&\cdots &[\,0\,]\\\vdots &\vdots &\ddots &\vdots \\[][\,0\,]&[\,0\,]&\cdots &{\hbox{ dens}}^{\,e}\end{bmatrix}} }}{\begin{matrix}1\\2\\\vdots \\n_{g}\end{matrix}}

(C.10)

And finally, the global density matrix is defined as:

\mathbf {\hbox{densglo}} ={\begin{bmatrix}{{\hbox{dens}}_{\,\mathbf {G} }}^{\,1}&[\,0\,]&\cdots &[\,0\,]\\[][\,0\,]&{{\hbox{dens}}_{\,\mathbf {G} }}^{\,2}&\cdots &[\,0\,]\\\vdots &\vdots &\ddots &\vdots \\[][\,0\,]&[\,0\,]&\cdots &{{\hbox{dens}}_{\,\mathbf {G} }}^{\,n_{el}}\end{bmatrix}}

(C.11)

After all this abstract formulation, the mass matrix can be simply assembled with no nested loops at all:

(C.12)

C.2 Generalisation to rotation matrices

During the development of the FEM method for rotating structures, it was found that new matrices emerged by the effect of rotation. These matrices have the form:

(C.13)

Where ${\textstyle \mathbf {A} }$ is a ${\textstyle n_{d}\times n_{d}}$ matrix accounting for the rotation effect, mostly Coriolis, centrifugal and tangential forces. The vectorised assembly of ${\textstyle \mathbf {P} }$ is almost straightforward, as most of the process is done exactly as in the case of the mass matrix. The only difference has to do with ${\textstyle \mathbf {A} }$ , which has to be globally defined:

{\mathbf {A} _{\mathbf {G} }}^{\,e}={\overset {(n_{d}\;n_{d}\;n_{g})}{\overbrace {\begin{bmatrix}\mathbf {A} ^{\,e}&[\,0\,]&\cdots &[\,0\,]\\[][\,0\,]&\mathbf {A} ^{\,e}&\cdots &[\,0\,]\\\vdots &\vdots &\ddots &\vdots \\[][\,0\,]&[\,0\,]&\cdots &\mathbf {A} ^{\,e}\end{bmatrix}} }}{\begin{matrix}1\\2\\\vdots \\n_{g}\end{matrix}}\;\;\;\;\;\;,\;\;\;\;\;\;\mathbf {A} _{\,glo}={\begin{bmatrix}{\mathbf {A} _{\mathbf {G} }}^{\,1}&[\,0\,]&\cdots &[\,0\,]\\[][\,0\,]&{\mathbf {A} _{\mathbf {G} }}^{\,2}&\cdots &[\,0\,]\\\vdots &\vdots &\ddots &\vdots \\[][\,0\,]&[\,0\,]&\cdots &{\mathbf {A} _{\mathbf {G} }}^{\,n_{el}}\end{bmatrix}}

(C.14)

Finally, the global matrix is assembled as:

(C.15)

D. Raw simulation results

In this chapter, results obtained from the developed simulations will be displayed. The discussion and analysis of the present results has already been addressed in the report (sections 6 and 7).

D.1 Two-dimensional cantilever beam

D.1.1 Natural modes and frequencies

Figure D.1: Natural modes associated to the eight lowest natural frequencies of a rectangular cantilever beam.

Table D.1. The twenty lowest natural frequencies of a rectangular cantilever beam.

D.1.2 Constant angular velocity

Case 1: Ω = 50 rad/s

Figure D.2: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.3: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Case 2: Ω = 100 rad/s

Figure D.4: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.5: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Case 3: Ω = 200 rad/s

Figure D.6: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.7: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Comparative figures

Figure D.8: Mode intensities for increasing values of

Figure D.9: Reaction torque for increasing values of

Figure D.10: Local transverse displacements at the tip of the beam for different values of

. Numeric results and its SVD approximation are compared.


(a)	(b)

(c)
Figure D.11: Power spectral density estimation of the sum of the right side vectors obtained for different rotation speeds: 50 (a), 100 (b) and 200 (c) rad/s.


(a)	(b)

(c)
Figure D.12: Local horizontal displacements at the tip of the beam for different values of $\Omega$ : 50 (a), 100 (b) and 200 (c) rad/s. Numeric results and its SVD approximation are compared.

Table D.2. Predominant frequencies and normalised singular values of the presented left side vectors.

Deformed shapes and stress distribution for Ω = 50 rad/s

Figure D.13: Deformed shape (

) and

distribution (colour map, in KPa) at 50 rad/s. Time from 0 to 14 ms.

D.1.3 Constant angular acceleration

Case 1: α = 5 rad/s²

Figure D.14: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.15: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Case 2: α = 50 rad/s²

Figure D.16: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.17: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Comparative figures

Figure D.18: Mode intensities for increasing values of

.

Figure D.19: Temporal evolution of vertical displacements at the tip of a cantilever beam under constant angular acceleration.

Figure D.20: Reaction torque for increasing values of

.

Table D.3. Predominant frequencies and normalised singular values of the presented left side vectors.

Deformed shapes and stress distribution for α = 50 rad/s²

Figure D.21: Deformed shape (

) and

distribution (colour map, in KPa) at

= 50 rad/s

and

= 50 rad/s. Time from 0 to 25 ms.

D.2 Three-dimensional cantilever beam

D.2.1 Natural modes and frequencies

Figure D.22: Natural modes associated to the eight lowest natural frequencies of a 3D cantilever beam.

Table. D.4 Fifty lowest natural frequencies of a 3D cantilever beam.

D.2.2 Constant angular velocity

Case 1: Ω = 25 rad/s

Figure D.23: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.24: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Case 2: Ω = 50 rad/s

Figure D.25: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.26: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Case 3: Ω = 100 rad/s

Figure D.27: Left side vectors (predominant modes) associated to the eight higher singular values.

Figure D.28: Right side vectors (predominant oscillations) associated to the eight higher singular values.

Campbell diagram

Figure D.29: Campbell diagram: comparison between modal prediction and SVD predominant frequencies. Subscripts

and

refer to pairs of equal modes that activate either on

or

.

Comparative figures

Figure D.30: Temporal evolution of

displacements at the tip of a cantilever beam under constant angular velocity.

Figure D.31: Temporal evolution of

displacements at the tip of a cantilever beam under constant angular velocity.

Figure D.32: Reaction torque for increasing values of

.

Deformed shapes and stress distribution for Ω = 50 rad/s

Figure D.33: Deformed shape (

) and

distribution (colour map, in KPa) at

= 50 rad/s. Time from 0 to 17 ms (up-down, left-right).

D.3 Starting of a rotor blade

Figure D.34: Temporal evolution of the angular speed and the local axial displacements at the tip.

Figure D.35: Temporal evolution of the angular acceleration and elastic reaction torque.

Figure D.36: Temporal evolution of the maximum values of stress.

Figure D.37: Right side vectors of the four most predominant oscillations.

Bibliography

[1] Rao, J. S. (2011) "History of rotating machinery dynamics". Springer 358

[2] Genta, Giancarlo. (2005) "Dynamics of Rotating Systems". Springer-Verlag New York

[3] Genta, Giancarlo, Silvagni, Mario. (2013) "On Centrifugal Softening in Finite Element Method Rotordynamics". Volume 81, Journal of Applied Mechanics

[4] Bentley, John R. (2007) "Aeolipile - Hero's Ball: A Replica of the World's First Rotating Steam Engine" http://modelengines.info/aeolipile/

[5] Volkov, Dale. (2016) "Heinkel He.178" http://www.airwar.ru/enc/xplane/he178.html

[6] Gao, Haifeng, Bai, Guangchen. (2015) "Reliability analysis on resonance for low-pressure compressor rotor blade". Advances in Mechanical Engineering http://journals.sagepub.com/doi/10.1177/1687814015578351

[7] Roylance, David "Beam Displacements". Department of Materials Science and Engineering, Massachusetts Institute of Technology

[8] González, José A., Abascal, Ramón, Park, K.C. (2014) "Partitioned analysis of flexible multibody systems using filtered linear finite element deformational modes". International Journal for Numeric Methods in Engineering, 102–128

[9] UPC ETSEIB "Numerical Factory - Finite Elements: a first course" https://numfactory.upc.edu/web/FiniteElements.html

[10] Tulloch Andrew. (2014) "Fast Randomized SVD - Facebook Research" https://research.fb.com/fast-randomized-svd/

[11] Inman, Daniel J. (2014) "Engineering Vibration". Pearson

[12] Belytschko, Ted, Lin, Wing Kam, Moran, Brian, Wiley, John. (2000) "Nonlinear finite elements for continua and structures". John Wiley & Sons

[13] Hernández, J.A., Oliver, J., Huespe, A.E. (2014) "Computational Homogenization of Inelastic Materials using Model Order Reduction". Monograph Nº M141, CIMNE

[14] D. Callister, William, Rethwisch, David. (1985) "Materials Science and Engineering: An Introduction". John Wiley and Sons

[15] "Airfoiltools - Airofil Database" http://airfoiltools.com

[16] V. Hutton, David. (2004) "Fundamentals of Finite Element Analysis". Mc Graw Hill

[17] Matsushita, Osami, Tanaka, Masato, Kanki, Hiroshi, Kobayashi, Masao, Keogh, Patrick. (2017) "Vibrations of Rotating Machinery - Basic Rotordynamics: Introduction to Practical Vibration Analysis", Volume 1 360

[18] Géradin, Michel, Cardona, Alberto. (2001) "Flexible Multibody Dynamics - A finite elements approach". John Wiley & Sons

[19] Japhet, Caroline, Cuvelier, Francois, Scarella, Gilles. (2015) "An efficient way to perform the assembly of finite element matrices - <hal-00931066v2>"

[20] Shabana, Ahmed A. (1997) "Flexible Multibody Dynamics: Review of Past and Recent Developments". Multibody System Dynamics 189-222

[21] GUWAHATI, NPTEL - IIT. (2014) "Theory of Rotor Dynamics" https://nptel.ac.in/courses/112103024/23

@@ Line 291: / Line 291: @@
 <div id='img-1.3'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
-|-style="padding:10px;"
+|-
-|[[Image:Draft_Samper_691820313-monograph-Campbell_example.png|450px|]]
+|style="padding:10px;"| [[Image:Draft_Samper_691820313-monograph-Campbell_example.png|450px|]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 1.3:''' Example of Campbell diagram <span id='citeF-6'></span>[[#cite-6|[6]]]
+| colspan="1" style="padding-bottom:10px;" | '''Figure 1.3:''' Example of Campbell diagram <span id='citeF-6'></span>[[#cite-6|[6]]]
 |}
@@ Line 300: / Line 300: @@
 To illustrate a little bit this section, Figure [[#img-1.3|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====
@@ Line 2,756: / Line 2,756: @@
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
 |-
-|style="padding:10px;"|[[Image:Draft_Samper_691820313-monograph-campbellnew.png|588px|Campbell diagram: comparison between modal prediction and SVD results.]]
+|[[Image:Draft_Samper_691820313-monograph-campbellnew.png|588px|Campbell diagram: comparison between modal prediction and SVD results.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" style="padding:10px;"| '''Figure 6.2:''' Campbell diagram: comparison between modal prediction and SVD results.
+| colspan="1" style="padding-bottom:10px;"| '''Figure 6.2:''' Campbell diagram: comparison between modal prediction and SVD results.
 |}
@@ Line 2,764: / Line 2,764: @@
 Recall that the obtained results shown in [[#D.1.2 Constant angular velocity|D.1.2]] are just '''predominant''' patterns of the dynamic response. Thus, the RSVs representing reigning vibrations are no longer pure harmonics, but mostly the result of a certain combination of natural frequencies. In most cases, they match almost perfectly to a sine or cosine function, and hence the predominant frequency is easy to determine. That is not the case of the RSVs associated to axial vibrations (natural modes 2 and 7 from Figure [[#D.1.1 Natural modes and frequencies|D.1.1]]), which are the result of more than one predominant frequency. For instance, the pure axial motion has a <math display="inline">0</math> Hz frequency associated to it, caused by the offset in displacements the centrifugal force induces. To visualise the frequencies involved in the vibration phenomena, Figure [[#img-D.11|D.11]] displays the '''power spectral density''' estimation of the sum of the right side vectors.
 ====Singular values====
@@ Line 4,292: / Line 4,291: @@
 <div id='img-D.8'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
 |style="padding:10px;"|[[Image:Draft_Samper_691820313-monograph-nsv_freq.png|400px|Mode intensities for increasing values of Ω]]
@@ Line 4,300: / Line 4,299: @@
 <div id='img-D.9'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
 |style="padding:10px;"| [[Image:Draft_Samper_691820313-monograph-React_comparison.png|400px|Reaction torque for increasing values of Ω]]
@@ Line 4,308: / Line 4,307: @@
 <div id='img-D.10'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 40%;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
-|style="padding:10px;" | [[Image:Draft_Samper_691820313-monograph-ydisp.png|Local transverse displacements at the tip of the beam for different values of Ω. Numeric results and its SVD approximation are compared.]]
+|style="padding:10px;" |[[File:Draft_Samper_691820313_8123_ydisp.png|800px]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" style="padding:10px;" | '''Figure D.10:''' Local transverse displacements at the tip of the beam for different values of <math>\Omega </math>. Numeric results and its SVD approximation are compared.
@@ Line 4,557: / Line 4,556: @@
 |style="padding:10px;"|[[Image:Draft_Samper_691820313-monograph-ydisp3d.png|420px|Temporal evolution of ϱ displacements at the tip of a cantilever beam under constant angular velocity.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" style="padding-bottom:10px;"| '''Figure D.30:''' Temporal evolution of <math>\varrho </math> displacements at the tip of a cantilever beam under constant angular velocity.
+| colspan="1" style="padding-bottom:10px;padding-left:10px;padding-rigth:10px;"| '''Figure D.30:''' Temporal evolution of <math>\varrho </math> displacements at the tip of a cantilever beam under constant angular velocity.
 |}
@@ Line 4,566: / Line 4,565: @@
 |style="padding:10px;"|[[Image:Draft_Samper_691820313-monograph-zdisp3d.png|480px|Temporal evolution of z displacements at the tip of a cantilever beam under constant angular velocity.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" style="padding-bottom:10px;"| '''Figure D.31:''' Temporal evolution of <math>z</math> displacements at the tip of a cantilever beam under constant angular velocity.
+| colspan="1" style="padding-bottom:10px;padding-left:10px;padding-rigth:10px;"| '''Figure D.31:''' Temporal evolution of <math>z</math> displacements at the tip of a cantilever beam under constant angular velocity.
 |}
@@ Line 4,577: / Line 4,576: @@
 | colspan="1" style="padding-bottom:10px;"| '''Figure D.32:''' Reaction torque for increasing values of <math>\Omega </math>.
 |}
 ====<span id='lb-D.2.2.6'></span>Deformed shapes and stress distribution for Ω = 50 rad/s ====
 <div id='img-D.33'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
 |style="padding:10px;"|[[Image:Draft_Samper_691820313-monograph-3dstress.png|383px|Deformed shape (×500) and σ<sub>ςς</sub> distribution (colour map, in KPa) at Ω₀ = 50 rad/s. Time from 0 to 17 ms (up-down, left-right). ]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" style="padding-bottom:10px;"| '''Figure D.33:''' Deformed shape (<math>\times 500</math>) and <math>\sigma _{\varsigma \varsigma }</math> distribution (colour map, in KPa) at <math>\Omega _0</math> = 50 rad/s. Time from 0 to 17 ms (up-down, left-right).
+| colspan="1" style="padding-bottom:10px;padding-left:10px;padding-rigth:10px;"| '''Figure D.33:''' Deformed shape (<math>\times 500</math>) and <math>\sigma _{\varsigma \varsigma }</math> distribution (colour map, in KPa) at <math>\Omega _0</math> = 50 rad/s. Time from 0 to 17 ms (up-down, left-right).
 |}
@@ Line 4,592: / Line 4,590: @@
 <div id='img-D.34'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
-|[[Image:Draft_Samper_691820313-monograph-xvel.png|600px|Temporal evolution of the angular speed and the local axial displacements at the tip.]]
+|[[Image:Draft_Samper_691820313-monograph-xvel.png|700px|Temporal evolution of the angular speed and the local axial displacements at the tip.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure D.34:''' Temporal evolution of the angular speed and the local axial displacements at the tip.
+| colspan="1" style="padding-bottom:10px;padding-top:10px;"| '''Figure D.34:''' Temporal evolution of the angular speed and the local axial displacements at the tip.
 |}
 <div id='img-D.35'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
 |[[Image:Draft_Samper_691820313-monograph-reactacc.png|600px|Temporal evolution of the angular acceleration and elastic reaction torque.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure D.35:''' Temporal evolution of the angular acceleration and elastic reaction torque.
+| colspan="1" style="padding-bottom:10px;padding-top:10px;"| '''Figure D.35:''' Temporal evolution of the angular acceleration and elastic reaction torque.
 |}
 <div id='img-D.36'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
 |[[Image:Draft_Samper_691820313-monograph-stress_blade.png|600px|Temporal evolution of the maximum values of stress.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure D.36:''' Temporal evolution of the maximum values of stress.
+| colspan="1" style="padding-bottom:10px;"| '''Figure D.36:''' Temporal evolution of the maximum values of stress.
 |}
 <div id='img-D.37'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
 |-
 |[[Image:Draft_Samper_691820313-monograph-rsvblade.png|600px|Right side vectors of the four most predominant oscillations.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure D.37:''' Right side vectors of the four most predominant oscillations.
+| colspan="1" style="padding-bottom:10px;padding-top:10px;"| '''Figure D.37:''' Right side vectors of the four most predominant oscillations.
 |}
-===BIBLIOGRAPHY===
+=Bibliography=
 <div id="cite-1"></div>
-'''[[#citeF-1|[1]]]''' Rao, J. S. (2011) "History of rotating machinery dynamics". Springer 358
+[[#citeF-1|[1]]] Rao, J. S. (2011) "History of rotating machinery dynamics". Springer 358
 <div id="cite-2"></div>
-'''[[#citeF-2|[2]]]''' Genta, Giancarlo. (2005) "Dynamics of Rotating Systems". Springer-Verlag New York
+[[#citeF-2|[2]]] Genta, Giancarlo. (2005) "Dynamics of Rotating Systems". Springer-Verlag New York
 <div id="cite-3"></div>
-'''[[#citeF-3|[3]]]''' Genta, Giancarlo and Silvagni, Mario. (2013) "On Centrifugal Softening in Finite Element Method Rotordynamics", Volume 81. Journal of Applied Mechanics
+[[#citeF-3|[3]]] Genta, Giancarlo, Silvagni, Mario. (2013) "On Centrifugal Softening in Finite Element Method Rotordynamics". Volume 81, Journal of Applied Mechanics
 <div id="cite-4"></div>
-'''[[#citeF-4|[4]]]''' Bentley, John R. (2007) "Aeolipile - Hero's Ball: A Replica of the World's First Rotating Steam Engine" http://modelengines.info/aeolipile/
+[[#citeF-4|[4]]] Bentley, John R. (2007) "Aeolipile - Hero's Ball: A Replica of the World's First Rotating Steam Engine" http://modelengines.info/aeolipile/
 <div id="cite-5"></div>
-'''[[#citeF-5|[5]]]''' Volkov, Dale. (2016) "Heinkel He.178" http://www.airwar.ru/enc/xplane/he178.html
+[[#citeF-5|[5]]] Volkov, Dale. (2016) "Heinkel He.178" http://www.airwar.ru/enc/xplane/he178.html
 <div id="cite-6"></div>
-'''[[#citeF-6|[6]]]''' Gao, Haifeng and Bai, Guangchen. (2015) "Reliability analysis on resonance for low-pressure compressor rotor blade". Advances in Mechanical Engineering http://journals.sagepub.com/doi/10.1177/1687814015578351
+[[#citeF-6|[6]]] Gao, Haifeng, Bai, Guangchen. (2015) "Reliability analysis on resonance for low-pressure compressor rotor blade". Advances in Mechanical Engineering http://journals.sagepub.com/doi/10.1177/1687814015578351
 <div id="cite-7"></div>
-'''[[#citeF-7|[7]]]''' Roylance, David "Beam Displacements. Department of Materials Science and Engineering, Massachusetts Institute of Technology"
+[[#citeF-7|[7]]] Roylance, David "Beam Displacements". Department of Materials Science and Engineering, Massachusetts Institute of Technology
 <div id="cite-8"></div>
-'''[[#citeF-8|[8]]]''' González, José A. and Abascal, Ramón and Park, K. C. (2014) "Partitioned analysis of flexible multibody systems using filtered linear finite element deformational modes". International Journal for Numeric Methods in Engineering 102&#8211;128
+[[#citeF-8|[8]]] González, José A., Abascal, Ramón,  Park, K.C. (2014) "Partitioned analysis of flexible multibody systems using filtered linear finite element deformational modes". International Journal for Numeric Methods in Engineering, 102&#8211;128
 <div id="cite-9"></div>
-'''[[#citeF-9|[9]]]''' UPC ETSEIB "Numerical Factory - Finite Elements: a first course" https://numfactory.upc.edu/web/FiniteElements.html
+[[#citeF-9|[9]]] UPC ETSEIB "Numerical Factory - Finite Elements: a first course" https://numfactory.upc.edu/web/FiniteElements.html
 <div id="cite-10"></div>
-'''[[#citeF-10|[10]]]''' Tulloch Andrew. (2014) "Fast Randomized SVD - Facebook Research" https://research.fb.com/fast-randomized-svd/
+[[#citeF-10|[10]]] Tulloch Andrew. (2014) "Fast Randomized SVD - Facebook Research" https://research.fb.com/fast-randomized-svd/
 <div id="cite-11"></div>
-'''[[#citeF-11|[11]]]''' Inman, Daniel J. (2014) "Engineering Vibration". Pearson
+[[#citeF-11|[11]]] Inman, Daniel J. (2014) "Engineering Vibration". Pearson
 <div id="cite-12"></div>
-'''[[#citeF-12|[12]]]''' Belytschko, Ted and Lin, Wing Kam and Moran, Brian and Wiley, John. (2000) "Nonlinear finite elements for continua and structures". John Wiley & Sons
+[[#citeF-12|[12]]] Belytschko, Ted, Lin, Wing Kam, Moran, Brian, Wiley, John. (2000) "Nonlinear finite elements for continua and structures". John Wiley & Sons
 <div id="cite-13"></div>
-'''[[#citeF-13|[13]]]''' Hernández, JA and Oliver, J and Huespe, AE. (2014) "Computational Homogenization of Inelastic Materials using Model Order Reduction". Cimne
+[[#citeF-13|[13]]] Hernández, J.A., Oliver, J., Huespe, A.E. (2014) "Computational Homogenization of Inelastic Materials using Model Order Reduction". Monograph Nº M141, CIMNE
 <div id="cite-14"></div>
-'''[[#citeF-14|[14]]]''' D. Callister, William and Rethwisch, David. (1985) "Materials Science and Engineering: An Introduction". John Wiley and Sons
+[[#citeF-14|[14]]] D. Callister, William, Rethwisch, David. (1985) "Materials Science and Engineering: An Introduction". John Wiley and Sons
 <div id="cite-15"></div>
-'''[[#citeF-15|[15]]]'''  "Airfoiltools - Airofil Database" http://airfoiltools.com
+[[#citeF-15|[15]]]  "Airfoiltools - Airofil Database" http://airfoiltools.com
 <div id="cite-16"></div>
-'''[[#citeF-16|[16]]]''' V. Hutton, David. (2004) "Fundamentals of Finite Element Analysis". Mc Graw Hill
+[[#citeF-16|[16]]] V. Hutton, David. (2004) "Fundamentals of Finite Element Analysis". Mc Graw Hill
 <div id="cite-17"></div>
-'''[17]''' Matsushita, Osami and Tanaka, Masato and Kanki, Hiroshi and Kobayashi, Masao and Keogh, Patrick. (2017) "Vibrations of Rotating Machinery - Basic Rotordynamics: Introduction to Practical Vibration Analysis", Volume 1 360
+[17] Matsushita, Osami, Tanaka, Masato, Kanki, Hiroshi, Kobayashi, Masao, Keogh, Patrick. (2017) "Vibrations of Rotating Machinery - Basic Rotordynamics: Introduction to Practical Vibration Analysis", Volume 1 360
 <div id="cite-18"></div>
-'''[18]''' Géradin, Michel and Cardona, Alberto. (2001) "Flexible Multibody Dynamics - A finite elements approach". John Wiley & Sons
+[18] Géradin, Michel, Cardona, Alberto. (2001) "Flexible Multibody Dynamics - A finite elements approach". John Wiley & Sons
 <div id="cite-19"></div>
-'''[19]''' Japhet, Caroline and Cuvelier, Francois and Scarella, Gilles. (2015) "An efficient way to perform the assembly of finite element matrices - <hal-00931066v2>"
+[19] Japhet, Caroline, Cuvelier, Francois, Scarella, Gilles. (2015) "An efficient way to perform the assembly of finite element matrices - <hal-00931066v2>"
 <div id="cite-20"></div>
-'''[20]''' Shabana, Ahmed A. (1997) "Flexible Multibody Dynamics: Review of Past and Recent Developments". Multibody System Dynamics 189-222
+[20] Shabana, Ahmed A. (1997) "Flexible Multibody Dynamics: Review of Past and Recent Developments". Multibody System Dynamics 189-222
 <div id="cite-21"></div>
-'''[21]''' GUWAHATI, NPTEL - IIT. (2014) "Theory of Rotor Dynamics" https://nptel.ac.in/courses/112103024/23
+[21] GUWAHATI, NPTEL - IIT. (2014) "Theory of Rotor Dynamics" https://nptel.ac.in/courses/112103024/23