2010:020
M A S T E R ' S T H E S I S
Evaluation of 3D Rotordynamics Capabilities within NX Nastran
Florian Thiery
Luleå University of Technology Master Thesis, Continuation Courses
Engineering mechanics
Department of Applied Physics and Mechanical Engineering
Division of Solid Mechanics
Acknowledgements
First of all, I would like to thank my supervisor Jan-Olov Aidanpää as well as my co-supervisor Marcus Sandberg for giving me the opportunity to carry out my Master’s thesis at the Solid Mechanics Department and Division of Functional Product Development in Luleå University of Technology. They have always been very kind and friendly when answering questions or having new ideas to carry through this project.
Then I would like to thank Christina Hamsch, from the International of- fice, for allowing me to stay one more semester to complete the Master’s de- gree. On the other side I would also like to thank the "Scolarité Office" from ENSEM (engineering school) who allowed students like me to stay longer in LTU to finish their Master’s thesis.
I would also like to express my sincere gratitude to my family, parents and grand-parents from France, who always supported me with my studies and agreed with the way I decided to pursue my education whatever the choices I made in the last 5 years.
Finally, I am very grateful to Amir Jourak, who helped me giving litera- ture and always shared happiness during coffee breaks; Matthias Wissemberg, friend of mine who made a lot to obtain the authorization to stay in Sweden;
I also send a warmly appreciation to Luisa Sa Vitorino, Martin Johansson without forgetting all the other students (exchange or Swedish) for their love and support during my whole studies.
Luleå, April 2010
Florian Thiery
Abstract
In today’s industry, companies within the jet engine industry are trying to make their product development more effective. One way is to use product models to enable upfront evaluation through simulation.
To model the dynamic behaviour of rotors, mainly beam and shell ele- ments are used in the common rotor-dynamics software. As the geometry becomes more complex, 3D rotor-dynamics can be used to solve vibration problems as well as the interactions with bearings and the non-rotating struc- ture. This 3D-modelling could also give a more realistic post-processing than the current software in the market and the dynamics of the whole engine could be assessed.
To try to solve these rotor-dynamics problems, NX7 was used, a commer- cial CAD/CAM/CAE software developed by Siemens PLM software. The solver used is NX Nastran Rotor-Dynamics. The main interest in this choice (with the 3D capabilities) is that both the CAD and CAE will be performed on the same platform and adjustments can be done faster if some problems in meshing or modelling occur.
The aim of this Master’s thesis is thus to develop a whole engine with the Computer Aided Design capabilities of NX7 and investigate the possibilities of doing rotor-dynamics analyses with NX Nastran Rotor Dynamics. First of all, simpler models were studied to validate the analyses and more com- plicated elements were added piece by piece until all the problems are solved to build the whole engine and analyse it.
The results presented in this Master’s thesis reveal several aspects of
rotordynamics analysis. First of all, a section is dedicated to the mode shapes
of the rotor for several common models. Then plots of Campbell diagram
are shown and compared with theoretical results. Finally, a mass unbalance
analysis is performed for all types of rotor and the displacement and rotation
response are plotted for some particular nodes. These displacement curves
are also compared with the Campbell diagrams to check if the simulation are run properly.
The main problem in this Master’s thesis is that the program is today too undeveloped for efficient 3D modelling of rotors. Indeed, it was found that rotating and stationary parts could not be coupled in the Campbell diagram.
Also a lot of commands and execution needs to be done manually and the possibilities of post-processing have to be developed to facilitate an effective evaluation of the results.
The future of this project will be to try to integrate optimization tools to these analyses, in order to have a better design and better dynamics results.
That is why the 3D rotordynamics simulation within NX are important: it
could enable to design a whole engine, study the dynamics of it and optimize
the whole system to come to a more efficient engine, all of this by using the
same software. This could result in time-savings and more efficient design of
engines in years to come.
Table of symbols
Symbol Description Unit
m Mass kg
d,C Viscous damping kg/s
K Stiffness N/m (kg/s
2)
Ω Rotational speed rad/s
F Force Newton (kg.m/s)
E Young’s modulus GPa
ρ Density kg/m
3ν Poisson’s ratio no unit
I Second moment of area m
4α Flexibility coefficient m/N
H Moment of Momentum N.s/m
M Moment N/m
J
pPolar moment of inertia kg.m
2J
dDiametric moment of inertia kg.m
2ξ Damping ratio no unit
x,y Position m
˙x, ˙y Velocity m/s
¨
x,¨ y Acceleration m/s
2e Eccentricity m
θ,ϕ Angle rad
˙θ, ˙ϕ Angular velocity rad/s
θ, ¨ ¨ ϕ Angular acceleration rad/s
2Nota Bene: The notation used here are general. They can be com-
pleted with a subscript, which could mean "rotation", "translation", "x-
direction",...
Contents
1 Introduction 10
2 Generalities 13
2.1 Rotor dynamics . . . 13
2.1.1 Linear Jeffcott rotor . . . 13
2.1.2 The gyroscopic effect . . . 15
2.1.3 Critical speeds and Campbell diagram . . . 21
2.1.4 Response to unbalance . . . 22
2.2 Finite element method . . . 23
3 The NX Nastran analysis 27 3.1 Theoretical foundation of rotor dynamics within Nastran . . . 27
3.1.1 Additional terms in the equations of motion . . . 27
3.1.2 Equation of motion for the fixed reference system . . . 28
3.1.3 The Steiner’s inertia terms in the analysis . . . 30
3.1.4 Real eigenvalue analysis . . . 31
3.1.5 The fixed system eigenvalue problem . . . 31
3.1.6 Solution interpretation . . . 32
3.2 Equation of motion for modal frequency response . . . 32
4 Modelling of the 2D-3D rotor in NX7 34 4.1 The middle span rotor (Jeffcott rotor) . . . 34
4.1.1 The geometric modelling of the rotor . . . 34
4.1.2 Meshing . . . 35
4.1.3 The boundary conditions . . . 42
4.2 The Overhung Rotor . . . 42
4.2.1 The geometric modelling of the rotor . . . 42
4.2.2 1D-Meshing . . . 43
5 Defining NX Nastran input for rotor dynamics 45 5.1 Modifications of the .dat file . . . 45
5.1.1 The file management section . . . 45
5.1.2 The executive control section . . . 46
5.1.3 The case control section . . . 46
5.1.4 The bulk data section . . . 47
5.2 Interpretation of rotor dynamics output . . . 48
5.2.1 The .f06 file . . . 49
5.2.2 The .op2 file . . . 49
5.2.3 The .csv file . . . 50
5.3 Procedure to perform an entire rotor dynamics simulation . . 51
6 Analysis of the results 53 6.1 The middle span rotor . . . 53
6.1.1 Quality of mesh and accuracy . . . 53
6.1.2 Campbell diagram . . . 54
6.1.3 Parameters influencing the results . . . 59
6.2 The Overhung rotor . . . 60
6.2.1 Mode shapes and Campbell diagram . . . 60
7 Solving the mass unbalance problem 63 7.1 Mass unbalance cards and modal frequency analysis . . . 63
7.1.1 Additional cards in the bulk data section . . . 63
7.1.2 Dynamic load . . . 64
7.2 Mass unbalance response of the Jeffcott rotor . . . 66
7.2.1 Synchronous analysis with internal damping . . . 66
7.2.2 Synchronous analysis with external damping . . . 70
7.2.3 Asynchronous analysis . . . 72
7.3 Mass unbalance response of the overhung rotor . . . 74
8 Modelling and analyses of a more realistic rotor 77 8.1 Modelling of the rotating part . . . 77
8.1.1 Meshing . . . 78
8.1.2 Boundary conditions . . . 81
8.2 Analysis of the modal complex eigenvalue solution . . . 82
8.2.1 Shaft without blades and no non-rotating structure . . 82
8.2.2 Shaft without blades and with the non-rotating structure 84 8.2.3 Shaft with blades and no non-rotating structure . . . . 85
8.3 Frequency response analysis of the rotor alone . . . 87
9 Conclusion and discussion 90
List of Figures
2.1 The mid-span rotor . . . 13
2.2 The unbalance disc on a shaft . . . 14
2.3 Simple view of the overhung rotor . . . 15
2.4 Forces and moments acting on the system in 3D . . . 16
2.5 Forces and moments acting on the system in 2D . . . 17
2.6 Geometric representation of z’ . . . 18
2.7 Example of Campbell diagram . . . 21
2.8 Response of an unbalance rotor as a function of the frequency 23 4.1 Revolution of the sketch of the rotor . . . 35
4.2 Type of solid elements . . . 37
4.3 Properties of the shell elements . . . 38
4.4 Connections of 1D elements at the end of the shaft (ref [2]) . 40 4.5 The "spider" rigid connection (in blue) . . . 41
4.6 The overhung rotor model . . . 43
4.7 General View of the 3D-rotor . . . 43
4.8 Cage surrounding the shaft with bearings . . . 44
5.1 Data entries of the rotor dynamic analysis . . . 49
5.2 Diagram of the entire rotor-dynamics procedure . . . 50
6.1 Translation and tilting mode of the mid-span rotor . . . 54
6.2 Example of a CHEXA(20) mesh and the eigenfrequencies . . . 55
6.3 Campbell diagram for the midspan rotor . . . 57
6.4 Damping values as function of the speed . . . 58
6.5 1st and 2nd modes shapes of the overhung rotor . . . 60
6.6 Campbell diagram for the overhung rotor . . . 61
7.1 Bulk data entries for the dynamic load . . . 65
7.2 Visualization of the dynamic load . . . 65
7.3 Displacement for the backward whirl . . . 67
7.4 Rotation for the backward whirl . . . 68
7.5 Theoretical displacement calculated with Matlab for the back- ward whirl . . . 69
7.6 One and a half degree of freedom system . . . 70
7.7 Displacement response for several damping ratios . . . 71
7.8 Displacement for the backward whirl of the asynchronous anal- ysis . . . 72
7.9 Rotation for the backward whirl of the asynchronous analysis . 73 7.10 Displacement of node 909 for the backward simulation . . . . 75
7.11 Displacement of node 909 for the forward simulation . . . 76
8.1 Rotor modelling with blades . . . 78
8.2 Real rotor to be modelled . . . 79
8.3 Meshing of the blades . . . 80
8.4 Connection of the bearings to the non-rotating structure . . . 81
8.5 Mode shapes of the rotor . . . 83
8.6 Unrealistic modes of the rotor . . . 84
8.7 First two modes of the rotor with the structure . . . 85
8.8 Mode shape of the blades . . . 86
8.9 Campbell diagram of the realistic rotor . . . 87
8.10 Backward and forward displacement of the realistic rotor . . . 88
Chapter 1 Introduction
The engineering design challenge presented by aerodynamic and hydrody- namic flows, design of blading, can sometimes cause the rotordynamics re- quirements of turbo-machinery (e.g. jet engines, gas turbines, compressors) design to be overlooked. It has happened that expensive turbomachines have been built and found to be incapable of producing their rated performance due to assumptions in the model.
Thus in designing, operating and troubleshooting turbo-machinery, rotor dynamics analysis can help accomplish the following objectives [3] :
1. Predict critical speeds (speeds at which vibration due to rotor imbal- ance is a maximum)
2. Determine design modification to change critical speeds
3. Predict amplitudes of synchronous vibration caused by rotor imbalance 4. Predict threshold speeds and vibration frequencies for dynamic insta-
bilities
5. Determine design modification to suppress dynamic instabilities
So basically, this project is considering three parts: designing, performing
rotor-dynamic analysis, upgrading the engine (optimizing), but in reality
more complex aspects to model have to be taken in account. The following diagram shows how they are related.
The part which will be developed in this project is the rotor dynamics analysis in 3D. Normally, these studies are done by using beam elements in the common rotor dynamics software (see Appendix A). Currently, only a few software are able to perform 3D rotor dynamics analyses. They usually are non-specific codes for rotor dynamics (i.e. they were not created to perform only rotordynamics analyses). These are:
• Ansys, the version 11 workbench and classic is capable of solving the rotordynamic equations (3-D/2-D and beam element)
• SAMCEF, a finite element based code (3-D/2-D and beam element)
• Nastran
Even if these software are not made to study the rotor-dynamics as a
priority, the most interesting part is that compared with the 1D beam element
solvers, all types of element can be used (beam, shell, solid elements) and
connected all together with springs and dampers to a big structure when the
whole designed model is very complicated. But still investigation on these
elements has to be done to be sure that the 3D rotordynamics analyses are
done properly with simple structures. Once this is secured, the 3D meshing
could be develop in a good way to build more accurate models.
Chapter 2 Generalities
2.1 Rotor dynamics
The aim of this chapter is to introduce the general concept of rotordynamics analyses. This is done by deriving the equations from a well-known model called the Jeffcott Rotor [4], also called the "De Laval rotor" in Europe.
2.1.1 Linear Jeffcott rotor
Figure 2.1: The mid-span rotor
The figure 2.1 shows a view of the general Jeffcott rotor. This model includes unbalance, viscous damping and whirl-instability (see section 2.3).
To understand the behaviour of the Jeffcott rotor, figure 2.2 shows the view of the disc with a mass unbalance. The disc is mounted on a shaft and is unbalanced, i.e. differs by the distance e from the geometric center (GC).
Figure 2.2: The unbalance disc on a shaft
The equations of motion can be derived from the previous figures. A position vector ~r is used from the origin O till the center of gravity of the disc. Thus:
~r = (x + e cos(Ωt))~i + (y + e sin(Ωt))~j (2.1)
~r = (¨ ¨ x − eΩ
2cos(Ωt))~i + (¨ y − eΩ
2sin(Ωt))~j (2.2)
The stiffness and damping of the shaft are called k and c respectively, and also m the mass of the disc. The Newton’s second law applied to the system gives:
m¨ ~r = −(kx + c ˙x)~i − (ky + c ˙y)~j (2.3)
Equation 2.2 into equation 2.3 gives the equations of motion:
m¨ x + c ˙x + kx = meΩ
2cos(Ωt) (2.4)
m¨ y + c ˙y + ky = meΩ
2sin(Ωt) (2.5) When ˙Φ = Ω, this gives the phenomena of synchronous whirl. This occurs when the angular velocity of the bended shape is equal to the angular velocity of the rotating speed. The whirl is a two dimensional orbit motion made by the geometric center. The orbit can have the same or the opposite direction as the spin direction of the rotor (respectively called forward whirl and backward whirl).
2.1.2 The gyroscopic effect
The gyroscopic effect happens when a spinning object is subjected to a rotational perturbation when a restoring moment appears. At high spin speeds or high moments of inertia, this phenomena becomes much significant.
This is analysed by determining the equations of motion of an overhung rotor.
Figure 2.3: Simple view of the overhung rotor
The stiffness of the beam can be obtained with the beam theory. The
following expressions can be found in an elementary case (cantilever beam
submitted to an end force and end moment).
For a beam submitted to a end force F, both displacement and phase are given by:
δ = F L
33EI = α
1F (2.6)
φ = F L
22EI = α
2F (2.7)
For a beam submitted to a end moment M, both displacement and phase are given by:
δ = ML
22EI = α
2M (2.8)
φ = ML
EI = α
3F (2.9)
The associated moments and forced can be derived from a free body diagram. A free body diagram is a pictorial representation of the forces acting on a free body, showing both contact and non-contact forces acting on this body (see figure 2.4).
Figure 2.4: Forces and moments acting on the system in 3D
If the deformations of the beam are studied in 2D with the equations (2.6)-(2.9), the flexibility matrix for the beam can be derived with the help of figure 2.5
Figure 2.5: Forces and moments acting on the system in 2D
Thus the following equations are obtained:
x = α
1F
x+ α
2M
y(2.10)
y = α
1F
y+ α
2(−M
x) (2.11)
−ϕ = α
2F
y+ α
3(−M
x) (2.12)
θ = α
2F
x+ α
3M
y(2.13)
In a generalized form, this gives:
q = [α]F =⇒ F = [α]
−1q = [K]q (2.14) The α matrix contains the flexibility coefficients like this:
[α] =
α
10 0 α
20 α
1−α
20 0 −α
2α
30
α
20 0 α
3
(2.15)
Thus the stiffness matrix ([K] = [α]
−1in equation 2.14) is given by:
[K] =
K
10 0 K
20 K
1−K
20 0 −K
2K
30
K
20 0 K
3
(2.16)
In this way the forces on the beam can be solved. The forces on the disc are the same but with opposite signs.
The gyroscopic effect appears when the beam bends. Indeed, the disc will rotate in a new direction Z’ which deviates from Z. The change in rotational direction results in a gyroscopic restoring moment of the disc. To describe and derive these moments, the moment of momentum H has to be determined.
The moment is equal to the time derivative of the momentum, i.e. ˙ H = M.
Thus it follows that:
H
z′= ΩJ
p(2.17)
Figure 2.6: Geometric representation of z’
From the geometry of figure 2.6, it follows that:
H
zz2′+ H
xz2 ′+ H
yz2 ′= H
z2′(2.18)
H
xz2 ′= H
zz2′tan θ (2.19)
H
yz2 ′= H
zz2′tan(−ϕ) (2.20) As the angles are assumed to be small, the equations 2.17, 2.19 and 2.20 into equation 2.18 give:
H
zz2′(1 + θ
2+ ϕ
2) ≃ H
z2′(2.21)
H
zz2′≃ ΩJ
p(2.22)
The equation 2.22 into equation 2.17 and 2.18 gives the components in x and y. According to figure 2.5, the moment of momentum from the rotations around ϕ and θ should be added. Thus one get:
H
x= J
dϕ + ΩJ ˙
pθ (2.23)
H
y= J
d˙θ − ΩJ
pϕ ˙ (2.24)
H
z= ΩJ
p(2.25)
The time derivative of the functions gives the moments caused by the gyroscopic effect:
M
Hx= J
dϕ + ΩJ ¨
p˙θ (2.26)
M
Hy= J
dθ − ΩJ ¨
pϕ ˙ (2.27)
M
Hz= 0 (2.28)
The equations 2.15, 2.26 and 2.27 give together with the result of rotating
unbalance the following equations of motion:
m¨ x = −K
1x − K
2θ + meΩ
2cos(Ωt) (2.29)
m¨ y = −K
1y + K
2ϕ + meΩ
2sin(Ωt) (2.30)
J
dϕ + ΩJ ¨
p˙θ = K
2y − K
3ϕ (2.31)
J
dθ − ΩJ ¨
pϕ = −K ˙
2x − K
3θ (2.32)
Equations 2.29 to 2.32 can be written into a matrix form:
m 0 0 0
0 m 0 0
0 0 J
d0 0 0 0 J
d
¨ x
¨ y
¨ ϕ θ ¨
+ Ω
0 0 0 0
0 0 0 0
0 0 0 J
p0 0 −J
p0
˙x
˙y
˙ ϕ
˙θ
+
K
10 0 K
20 K
1−K
20 0 −K
2K
30
K
20 0 K
3
x y ϕ θ
= meΩ
2
cos(Ωt) sin(Ωt)
0 0
(2.33)
In a shorter way, it can also be written:
[M]¨ r + Ω[G] ˙r + [K]r = [F ] (2.34)
The equation can be used to simulate the dynamics of the system, for ex-
ample by using Matlab. The property of the gyroscopic matrix is dependent
of Ω, that means the eigen-frequencies will depend on the rotating speed of
the rotor. This is a specific result of the rotating systems.
2.1.3 Critical speeds and Campbell diagram
Every rotor-bearing system has a number of discrete natural frequencies of lateral vibration. When one of the natural frequencies is excited by rotor im- balance rotating at shaft speed, the shaft speed which coincides with natural frequency is called a critical speed.
In mathematical terms, the natural frequencies are called eigenvalues and the mode shapes eigenvectors. The mode shape describes the expected curvature (or displacement) of the rotor at a particular mode. Theoretically, a distributed mass-elastic system has an infinite number of eigenvalues and associated eigenvectors. But in practice, only the lowest critical speeds and associated whirl modes are excited in the speed range of a typical turbo- machine.
Figure 2.7: Example of Campbell diagram
As the natural frequencies vary with the spin speed (equation 2.34), the best way to find these critical speeds is to plot the Campbell diagram [5].
It consists of plotting the natural frequencies as a function of the rotational speed. Theoretically, these damped natural frequencies are the imaginary part of the eigenvalues. The critical speeds are found where the natural frequencies coincide with the speed spin (ω = Ω). In figure 2.7, frequencies that increase (blue line) is called forward whirl whereas the one that decrease (red line) is named backward whirl. The straight red line shows a linear whirl, which means that the eigen-frequency is not depending on the rotational speed. The critical speeds can be observed when the black line (called 1P line) crosses the other curves (red or blue).
2.1.4 Response to unbalance
As shown in equation 2.33, a sinusoidal force is applied on the system. This is due to the unbalanced mass, which is one of the main problems in rotor- dynamic systems. For a 1-dimensional system, the equation to solve is:
m¨ x + c ˙x + kx = m
ueΩ
2sin(Ωt) (2.35) The excitation is due to an applied force F
0sin(Ωt) where F
0= m
ueΩ
2. This is the centripetal force of the mass m
utowards the center of rotation.
If the rotor is balanced, no force exists, but this is very seldom the case. The amplitude of the excitation depends on both the out-of-balance mass, and its effective distance from the axis of rotation, as well as on the rotational speed.
The response will be given by x = X
0sin(Ωt − φ) where:
X
0= m
ueΩ
2/k
p(1 − r
2)
2+ (2ξr
2) (2.36) tan φ = 2ξr
1 − r (2.37)
Here ω is the natural frequency ω = r k
m , ξ = c 2 √
km and r = Ω
ω .
In a non-dimensional form, it is easier to see how in general the response depends on the frequency ration r and on the damping ratio ξ. Therefore:
X
0e
m m
u= r
2p(1 − r
2)
2+ (2ξr
2) (2.38)
Then the figure 2.8 is obtained by plotting X
0e
m m
uagainst r for several damping ratios.
Figure 2.8: Response of an unbalance rotor as a function of the frequency
2.2 Finite element method
The finite element method is a very general method to handle difficult prob-
lems that cannot be solved analytically. This method is widely used in many
fields (thermodynamics, electrical engineering, ...) and one that is very suit- able for is in structural mechanics. FEM is based on the subdivision of the structure into finite elements.
Element formulation
Many different element formulations have been developed depending on their shape and characteristics: beam elements, shell elements, plate elements, solid elements and many others. Each element is the model of a small de- formable solid. Inside each element the displacement u(x, y, z) of the point with the coordinates x,y,z is approximated by combination of the shape func- tions N. FEM is usually developed using matrix notation to obtain formulas that can easily be transferred into computer code. The displacements are usually written in three dimensions:
u(x, y, z, t) = N(x, y, z)q(t) (2.39)
The displacements can contain more dimensions if rotations of the el- ements are considered. In this case, q is a vector with the n generalized coordinates, and N the matrix containing the shape functions. There are as many rows in N as in u and as many columns as the number n of degrees of freedom. For a two dimensional element, the equations can be written:
u
x(x, y, t) u
y(x, y, t)
!
= N(x, y) 0
0 N(x, y)
! q
x(t) q
y(t)
!
(2.40)
The shape functions can be chosen arbitrarily, but they have to satisfy several conditions such as:
• be continuous and derivable up to the required order, which depends on the element type
• be able to describe rigid body of the element leading to vanishing elastic
potential energy
• lead to constant strain field when the overall deformation of the element dictates so
• lead to deflected shape of each element that matches the shape of the surrounding elements
The strains in each element can be expressed as functions of the deriva- tives of the displacements with respect to space coordinates; this makes pos- sible to write the strains as:
ε(x, y, t) = B(x, y)q(t) (2.41) B is a matrix containing the derivatives of the shape functions. The stress can be expressed as:
σ(x, y, t) = Eε = E(x, y)B(x, y)q(t) (2.42) E is the stiffness matrix of the material. Since only isotropic material is used, E is a scalar. The potential energy can then be expressed as:
U = 1 2
Z
V
ε
TσdV = 1 2 q
T(
Z
V
B
TEBdV )q (2.43)
The integral in equation 2.43 is the stiffness matrix of the element:
K = Z
V
B
TEBdV (2.44)
The mass matrix is derived from the expression of kinetic energy of the element and the shape functions:
T = 1 2 ˙q
T(
Z
V
ρN
TNdV ) ˙q (2.45)
M = Z
V
ρN
TNdV (2.46)
The form of the matrix B is dependent of the shape functions N and
can hence be in different dimensions and considered in different degrees of
freedom. This gives different dimensions of the stiffness matrix depending of the assumed shape functions. Then these two matrices are introduced into Newton’s second law:
M ¨ q + Kq = F (t) (2.47)
The Fung and Tong book [6] can be investigated to have more details
about the finite element method (shape functions for solid elements, assembly
process, ...).
Chapter 3
The NX Nastran analysis
This chapter has similarities with Chapter 2 as it presents theoretical foun- dation of rotor-dynamics.
3.1 Theoretical foundation of rotor dynamics within Nastran
A said in Chapter 1, for a rotating structure, additional terms occur in the equations of motion depending on the chosen analysis system. The fixed system will mainly be used, and the rotor is rotating about the z-axis (default axis). The equations of motion are described in this system. One has to set
"GRID" points that belongs to the rotor to distinguish between the rotating and fixed parts of the structure.
3.1.1 Additional terms in the equations of motion
• Coriolis forces and gyroscopic moments
Rotor dynamics analyses in NX Nastran include both Coriolis forces
and gyroscopic effects. In the rotating analysis system, the Coriolis
forces of the mass points are included. In the fixed analysis system,
the gyroscopic moments due to nodal rotations are included.
• Damping
The damping in a rotor system is divided into two parts: the internal damping acting on the rotating part of the structure, and the external damping acting on the fixed part of the structure and in the bearings.
The rotor can become unstable at speeds above the critical speed as the internal damping has a destabilizing effect. On the contrary, the exter- nal damping has a stabilizing effect (more information about damping is given in Appendix C).
3.1.2 Equation of motion for the fixed reference system
The generalized equation of motion is as follows:
[ ¯ M ]¨ q + (Ω[ ¯ C] + [ ¯ D
I+ ¯ D
A]) ˙q + ([ ¯ K] + Ω[ ¯ D
B])q = 0 (3.1)
The generalized matrices are given by:
• [ ¯ M ] = [Φ]
T[M][Φ] = [I] the generalized mass matrix
• [ ¯ C] = [Φ]
T[C][Φ] the generalized antisymmetric gyroscopic matrix
• [ ¯ D
I] = [Φ]
T[D
I][Φ] the generalized internal viscous damping matrix
• [ ¯ D
A] = [Φ]
T[D
A][Φ] the generalized external viscous damping matrix
• [ ¯ K] = [Φ]
T[K][Φ] = diag[ω
02] the generalized elastic stiffness matrix
• [ ¯ D
B] = [Φ]
T[D
B][Φ] the generalized anti symmetrical internal damping matrix
For a rotating mass point, the terms are as it will follow.
The mass matrix is the same as for non-rotating structures. A lumped
mass approach is used here:
[M] =
m
m m
J
xJ
yJ
z
(3.2)
In the gyroscopic matrix, only the polar moment of inertia appears. With- out any polar moments of inertia,there is no rotational effects in the fixed system.
[C] =
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 J
z0 0 0 0 −J
z0 0
0 0 0 0 0 0
(3.3)
The internal rotor damping matrix is given by:
[D
I] =
d
I,xd
I,yd
I,zd
I,Rxd
I,Ryd
I,Rz
(3.4)
The external damping acting in the non-rotating bearings is as follows:
[D
A] =
d
A,xd
A,y0 0
d
A,Rxd
A,Ry
(3.5)
Because the displacements in the rotating part act as a velocity in the fixed part, an antisymmetric matrix appears in the stiffness term :
[D
B] =
0 d
I,T0
−d
I,T0 0
0 0 0
0 d
I,R0
−d
I,R0 0
0 0 0
(3.6)
In the previous matrices, the following notations have been used:
m mass
J
xinertia about the x-axis J
yinertia about the y-axis J
zinertia about the z-axis d
Iviscous damping of rotor d
Aviscous damping of rotor
The subscript T refers to "translation" whereas the subscript R indicates
"rotation". Moreover, in equation 3.6, the damping is assumed to be equal in the x- and y-direction.
3.1.3 The Steiner’s inertia terms in the analysis
The ZSTEIN term is an option which can be found on the ROTORD card (the rotordynamic entries are defined in this card) in the data file. This allows to include the Steiner’s in the analysis. It is important to include this term in the fixed reference system. In this case, the polar moment of inertia is calculated as :
J
z= J
p= Σm(dx
2+ dy
2) (3.7) WARNING: to use the ZSTEIN option, the nodes must have rotations.
In a model made of solid elements, this is obtained by adding a layer of
shell elements around the solid elements. One should also check if the nodal rotations are not constrained by any optional cards.
3.1.4 Real eigenvalue analysis
The first step in a rotor dynamics analysis is to perform a real eigenvalue analysis :
(−ω
02[M] + [K])ϕ = 0 (3.8)
The eigenvectors ϕ are collected into the modal matrix Φ. The dis- placement vectors can be described by a linear combination of the modes:
u = [Φ]q. This is mathematically true if all modes are used. It is reasonable when using a sufficient number of modes. The selected modes used for the rotor dynamic analysis must be chosen according to these criteria:
• The frequency range of the real modes must be well above the frequency range of interest in the rotor dynamic analysis
• The real modes should be able to describe the rotor motion. The modal displacement in x and y has to be included
• The modes must be able to represent the generalized forces
3.1.5 The fixed system eigenvalue problem
For harmonic motion q(t) = q exp(λt), the following problem is solved in the fixed reference system:
(λ
2[ ¯ M ] + λ(Ω[ ¯ C] + [ ¯ D
I+ ¯ D
A]) + ([ ¯ K] + Ω[ ¯ D
B]))q = 0 (3.9)
A loop is made over the selected rotor speeds and the results are stored
for post processing.
Synchronous analysis
The eigenfrequency is equal to the rotor speed for the points crossing the 1P-line, i.e. ω = Ω. If the damping is neglected, λ = jω. To obtain the resonance points, the following eigenvalue equation is solved:
(−Ω
2([ ¯ M ] + j[ ¯ C]) + [ ¯ K])q = 0 (3.10) The imaginary parts of the solutions are the critical rotor speeds. For the modes that do not cross the 1P line, the imaginary part is zero.
3.1.6 Solution interpretation
The solutions at each rotor speed are the complex conjugate pairs of eigen- values λ = δ ± jω. The real part is a measure of the amplitude amplification.
Having positive values leads to an increase in amplitude with time and so the mode is unstable. The system is considered as stable when the real part is negative. If the damping is assumed low, it can be defined as:
ξ = δ
ω (3.11)
In the complex modal analysis within NX Nastran, the output is the g-damping defined as g = 2ξ.
The imaginary part of the eigenvalue is the damped natural frequency of the solution. The eigenfrequency is thus f = ω
2π . The units that can be chosen in Nastran will be shown in Chapter 4.
3.2 Equation of motion for modal frequency re- sponse
Fixed reference system
The governing equation of the frequency response in modal space with rotor
dynamics terms, in the fixed reference system considering the load to be
independent of the speed of the rotation is:
(−ω
k2[ ¯ M] + jω
k(Ω[ ¯ C] + [ ¯ D
I+ ¯ D
A]) + ([ ¯ K] + Ω[ ¯ D
B]))u(ω
k) = p(ω
k) for k = 1, 2, ..., m
Here, m denotes the number of excitation frequencies of dynamic load.
This is called an asynchronous solution and applicable to cases like gravity loads. The load could still have frequency dependence as shown by the m discrete excitation frequencies defined on the FREQ card (see Chapter 7).
The rotor speed is constant and the asynchronous analysis is working along a vertical line in the Campbell diagram.
In the case the load is dependent on the speed of rotation (called syn- chronous analysis) is found putting ω = Ω. The governing equation is as follows:
(−Ω
2k([ ¯ M ] − j[ ¯ C]) + jΩ
k([ ¯ D
I+ ¯ D
A] − j[ ¯ D
B]) + [ ¯ K])u(Ω
k) = p(Ω
k) for k = 1, 2, ..., n
Here n denotes the number of Ω
krotation speeds at which the analysis
is executed. Such is the case for example with centrifugal loads due to mass
unbalance. In this case the analysis is done along the 1P excitation line.
Chapter 4
Modelling of the 2D-3D rotor in NX7
In this section, two rotors will be modelled in order to validate the models created in NX 7. The results will be then compared with the theory if possible. The rotor which will be analysed are a mid-span rotor and an overhung rotor. The choice of both rotor has been done to try to model a range of rotors as different as possible.
4.1 The middle span rotor (Jeffcott rotor)
4.1.1 The geometric modelling of the rotor
The parameters of the mid-span rotor are L = 1 m, ρ = 7800 kg/m
3, E = 220 GPa, D = 1 m, d = 0.05 m and t = 0.1 m.
First of all, the model needs to be created under NX7. One have to open a new part and start to "sketch" the model. Once the sketch is completed, the solid body should be created by revolving the sketch around the z-axis.
The result of these operations gives the rotor figure 4.1.
Nonetheless, the design has to be studied before meshing. Indeed the
rotor can be modelled in two ways: either the rotor is made in one single
Figure 4.1: Revolution of the sketch of the rotor
part, or the shaft and the disc are made separately and glued together. The choice of the model has an influence on the mesh. For reasons of accuracy, the 2-part rotor will be used in the rest of the modelling.
4.1.2 Meshing
The mesh of the rotor is done in the Advanced Simulation application of NX7. To be more accurate, a FEM file should be created and it is in this section that the entire mesh will be completed. This section will be split in categories (1D, 2D and 3D meshes), but in the order used to create the model, i.e. the 3D mesh first, then the 2D mesh and the 1D to finish.
4.1.2.a The 3D mesh
Three important points have to be considered when meshing the rotor:
• The element type
• The element size
• For a 2-part rotor, the mesh mating condition between the contact surfaces
The element type
Fundamentally, there are two types of 3D mesh which could be used. There are referenced as CTETRA and CHEXA. The first one is a tetrahedron whereas the second one is a 6-face element. Moreover in each element one can choose either an element with nodes just in the corners or an element with nodes in the corners and in the middle of the edges. For the CTRETA element, it can be CTETRA(4) or CTETRA(10), and for CHEXA one can pick CHEXA(8) or CHEXA(20). The number in parenthesis indicates the number of nodes for each element. It is important to choose elements with nodes on the edges because the accuracy will be much better using them.
WARNING: The elements with the greater number of nodes should always be used. Indeed, the accuracy of the elements with only nodes in the corner is really poor when simulating. For instance, the strains are constant in a 4-node tetrahedron, which is the simplest solid element and which is not very accurate [6].
In addition to this, a CPENTA element can sometimes be found in the data file. This is an automatic creation by the software in some parts of the rotor when the CHEXA-mesh is used.
After completing the 3D mesh, it is important to choose a material in the properties of the 3D elements, e.g. steel, or create his own material (see section 4.1.2.d), otherwise the simulation will not be launched.
The element size
The size of the element chosen throughout the body has an effect on the
overall results. If one select a too large size, this can result in a rotor which
is farther from the initial model. Indeed the shaft of the rotor can look more
like a polygonal object than a cylindrical one. To make a first test, the best
Figure 4.2: Type of solid elements
way is to try an action called "Automatic Element Size" in Nx. This is a good approach to start with, but then it can be modified manually to find the best results. Then one can choose smaller elements through the shaft and the span, but the effect will be that the fastness of the simulation will drop down as the size of the element decreases. Therefore a compromise should be found to have both good results and a rather good quickness of the simulation.
The mesh mating condition
As the rotor can be divided in several parts, it is important that the meshes should be well-connected, the nodes being the same between the contact sur- faces. This is done by using a command in the fem file called "Mesh Mating Condition". If this is not done, the simulation will result in displacements of both the shaft and the span totally disconnected, being independent from each other.
4.1.2.b The 2D mesh
For this example, the 2D mesh will be used just for one goal. Indeed, to have
the analysis in the fixed reference system, it is necessary to put the ZSTEIN
entry to YES in the rotor bulk section (cf Chapter 5) and that the nodes have
rotations. In a model with solid elements, this can be done by adding a layer
of shell elements around all the solid elements.These shell elements have a
negligible effect on the model’s overall stiffness. When using this option, one must ensure that the nodal rotations of the shell elements are not constrained by the AUTOSPC option.
To insert all these shell elements in NX in the advanced simulation win- dow, the path to follow is: Insert->Element->Surface Coat. Then all the solid elements should be selected. There is no need to choose the type of element (i.e. CQUAD8 or CTRIA6) as the creation will be automatic.
Once this is done, the properties of the shell elements has to be completed.
The two most important ones are Material 1 and Default Thickness. In this model, the material will be steel and the shell elements will be 1 mm thick.
Figure 4.3: Properties of the shell elements
4.1.2.c The 1D mesh
The 1D mesh elements will be found at both ends of the rotor. All these 1D
elements will be connected to each other to the central point of the rotor.
But as the software does not allow to connect for example one spring from one node to the same, it is useful to create 3 nodes in the middle of each end surface.
As the connection between the 3D elements and 1D may cause some problems, it is important to insert a rigid connection between the bearings and the end surface. Then the bearings will be attached to this connection and a end constrained point. To sum up, there are 3 types of 1D connections in the model:
• The rigid connection (RBE2)
• The bearings (CELAS1)
• The damping of the bearings (CDAMP1)
The rigid connections
There are two types of rigid connections inserted into the model: one between all the points of the end surface and the middle node of this same surface and another connection from this middle node to the bearing node (both nodes have the same coordinates).
The first connection is called "spider connection". One have to go in Insert->Mesh->1D Connection, to choose Node to node in the connection type box, to select RBE2 in the type of element. Then the central node has to be selected as a source and all the other nodes of the surface should be selected as a target. Plus the rotational degrees of freedom have to be let free to allow rotations in X,Y and Z. Doing that, there will be no problem of connection while simulating.
The second connection is made in the same way as the previous one,
excepted that the source is the same but the target node is one at the same
coordinates which has been created before. Plus all degrees of freedom need
to be constrained in this RBE2 connection.
Figure 4.4: Connections of 1D elements at the end of the shaft (ref [2])
The most important thing to remember is that one should be careful in
which sense to connect these nodes, otherwise over-constrained nodes may
appear when putting the end constraints.
Figure 4.5: The "spider" rigid connection (in blue)
The bearings and damping of bearings
The bearings and damping of bearings need to be inserted between the second node and the third node created in the middle of the end surface. As the RBE connection, it can be found in Insert->Mesh->1D Connection by changing the type of element. In this model, two bearings will be created as well as damping of bearings in the X and Y-direction.
For the viscous damping (external), a small value has been chosen in order to avoid numerical problems as a start. The value has been fixed to 10
7N.m for the stiffness.
4.1.2.d Material
There is two choices concerning the choice of the material: either one pre-
defined material can be selected (such as Isotropic Steel ), or a new material
can be created. It is better to create an own material as the inserted values
would be then the real ones. It is also in this case that one can put some
value for the internal damping, as the predefined material do not have it (and
the values of the predefined material cannot be changed).
4.1.3 The boundary conditions
As shown on figure 4.4, both end nodes will be totally constrained (i.e. the six degrees of freedom). To do that, one have to go in User Defined Constraint in the simulation file. Then the nodes must be selected one by one and all degrees of freedom have to be fixed.
Then the end nodes of the rotor should be constrained as the beam is simply supported. That means one have to proceed as before, but only the displacements in the x,y and z-directions are fixed as well as the rotation around the z-axis.
Once all of this done, the first step is to analyse the eigen-frequencies of the rotor to see if the mesh is satisfying. Then the rotor dynamics entries can be inserted into the ".dat" file.
4.2 The Overhung Rotor
The aim of studying this rotor is that unlike the mid-span rotor, the eigen- frencies of the first modes will not be the same as the rotor speed increases, because of the gyroscopic effects. Due to that, all types of simple rotors can be studied whatever their geometry. Plus as the boundary conditions are changed, one have to find another way to put the bearings in this model without affecting the overall results.
4.2.1 The geometric modelling of the rotor
The parameters of the mid-span rotor are L = 1 m, ρ = 7800 kg/m
3, E = 220 GPa, D = 1 m, d = 0.05 m and t = 0.1 m (see figure 4.6).
As the mid-span rotor, the overhung rotor will me modelled with two parts, i.e. the shaft and the disc. Then the control of the size of the mesh is much easier.
Once the sketch is done regarding the dimensions, one have to revolve it
to finally obtain the 3D-rotor (figure 4.7)
Figure 4.6: The overhung rotor model
4.2.2 1D-Meshing
As the 3D-mesh and 2D-mesh will be the same as in section 4.1.2, the 1D mesh will be emphasized as the boundary conditions differ from the mid-span rotor.
Indeed, it is impossible to put the bearing as in the mid-span rotor because it will not be taken in the analysis proceeding in this way. Thus the idea is to create a small cage around the end of the rotor in order to be able to connect the bearings from this "cage" to the rotor.
Figure 4.7: General View of the 3D-rotor
But still some problems can appear by working like this. How many nodes have to be created to be sure that the boundary condition are respected? How many bearings should be created and to which rotor nodes should they be connected? What is the value of each bearing as several are connected to each other?
The answers to these questions will be analysed and discussed with the results of the simulation. But to start with something, the idea is to create three bearings on each side of the small outside cage and connect them to two center points of the rotor (one for each side).
Figure 4.8: Cage surrounding the shaft with bearings
If the results are satisfying, this model will be considered as a good one, otherwise some elements will be changed. Another important thing is that small damping will be introduced first to avoid numerical issues, trying to be as more time-saving as possible.
Concerning the 2D and 3D mesh, it is exactly the same for the Jeffcott
rotor.
Chapter 5
Defining NX Nastran input for rotor dynamics
The aim of this chapter is to describe the main commands to properly analyse the rotor and structure and show how to process the results. Indeed, as described in chapter 3, the solution used is the analysis of modes also known as SEMODES-103. But to run the rotor dynamic analysis, one need to run solution 110 (complex modal) which is not supported by NX7. Thus once the solution has been launched, one have to Write, Edit and Solve in the Solution Attribute Box and save the text file as a .dat file. Then modifications can be done manually by opening the .dat file and adding the bulk data entries of the rotor as well as changing the file management section, the executive control section and the case control section.
5.1 Modifications of the .dat file
5.1.1 The file management section
In the file management input file, one should assign statements to have results
files which will be used for post-processing. These files generated by the
software are called GPF and CSV. To get it, the followings statements need
to be written in this section:
assign output4=’name_of_file.gpf’, unit=22, form=formatted
assign output4=’name_of_file.csv’, unit=25, form=formatted
