Discretization techniques

According to Genta (2005), three dierent classes of the discretization

tech-niques can be distinguished. The rst class of the discretization techniques

comprisesassumed mode methods. These methods assumethat the deected

shapeofacontinuouselasticbodyisalinearcombinationof

M

^{arbitrarymode}

shapefunctions

φ _k ( x )

^(for

k ∈ [1, M ]

^)ofthespacecoordinates

x

^{. Asan}

exam-ple, thedeection

y _j

^{of aspecic pointof thebody (in theundeformed body}

thepoint

x _j ∈ V

^,wherethe

V

^{isthespatialvolumeoccupiedbythebody)at}

time

t

^{canbeexpressedintermsoftheassumed}

φ _k (x j )

^{andmodalamplitudes}

η k (t)

^as

where

η k

^{can beconsideredasthemodal}coordinates. Furthermore,weassume that the displacement of the structure is measured at the

P

^{discrete points}

(nodes)

x = x j

^(for

j ∈ [1, P ]

^{)andthereforethe}measurementvectoreld is

Inthematrixnotation,Eq.(3.13)canbeexpressedas

y = Φη,

^(3.14)

Fromtheobtaineddisplacementeld,theexpressionsofthekinetic,potential,

and dissipationenergies canbe obtained. The totalkineticenergy ofthe

sys-tem (withthe positive anddenite massmatrix, that is, for instance, forreal

symmetric matrix

M

^{and for all non-zero real vectors}

η ˙

^,

η ˙ ^T M η ˙ > 0

^{) is}

ex-pressed in termsof the generalizedcoordinates

η

^and generalizedvelocities

η ˙

as

W T = W T (η, η) ˙ .

^{Thepotentialenergyis expressedintermsofthe}

general-izedcoordinatesonly

W V = W V (η)

^andthedissipationenergydependsonthe generalized velocities only

W P = W P ( ˙ η)

^{. In general, for our purposes, these}

expressionscanbeexpressedas

W T = 1

2 η ˙ ^T M η, W ˙ V = 1

2 η ^T K M η, W P = 1

2 η ˙ ^T C M η, ˙

^(3.16)

where

K M

^and

C M

^{arethegeneralstinessmatrixandgeneraldampingmatrix}

that consists of a skew-symmetric gyroscopic matrix

G

⁽

G ^T = − G

^{) and a}

symmetric damping matrix

D M

^{. The}

M

^,

K M

^and

C M

^{are of order}

M

^and

they depend on the inertial and elastic properties of the system. Now, the

equations of motion can becomputed using Lagrange'sequation. Lagrange's

equationforaparticularDOF(associatedwiththe

η k

^{)canbewritten}

d

isthemodalforceequaltothesumofthediscreteforcecomponents

f j

^applied

at each node (

x j

^{). Theequations of motion, in modal}coordinates

η

^{, derived}

from Lagrange's equations for the rotational speed

Ω

^{can be written in the}

matrixform

M η ¨ + (D M + ΩG) ˙ η + K M η = f .

^(3.19)

The second group of the discretization techniques are lumped-parameter

methods. Inthesemethods,thephysicalstructureis dividedintoanumberof

rigid bodies connectedbymasslesseldsthat possesstheelastic anddamping

properties. Here,the generalizedcoordinatesaredened bythedisplacements

of therigidbodies. Furthermore,using thesameapplication ofNewton's

Sec-ond Lawasforthe 1-DOF model, wecancompose aM-DOFmodel. The

M

and

f

^{are written directlyusing the}engineeringjudgment. However, for com-plex structures, obtainingthe

K M

^{matrix is dicult. Therefore, the}

lumped-parametermethodsoftenresorttonite elementmethods(FEMs) tocompute

thestinessmatrix.

The third most versatile class of the discretization techniques are FEMs

(alsoassumed modes areusually derived usingFEM), which arebasedon the

subdivision ofthestructure into niteelements. Each elementhasits ownset

of DOF associated with the discrete nodes. Inside each element, similarly as

in the assumed mode method, the displacement

w(x, t)

^is approximated by the usually linearcombination of the shapefunctions (Eq. 3.13),where

η

^are

nowthegeneralizedcoordinatesoftheparticularelement. Theshapefunctions

φ _k

^{are arbitrary but must satisfy a number of} conditions, particularly to be

continuousandderivable,andtheyhavetomatchtheshapeoftheneighboring

elements. Theneighboring elementsshare the nodesplaced at theside of the

elements. Foreachelement,theequationsofmotioncanbewrittenanalogically

asfor theassumedmodesmethod. Next,thecomplete setof equationshasto

be assembled. In thecase of therotormodeling, thedisk and beamelements

canbeused; iftheirreferenceframes coincidewith theglobalreferenceframe,

nocoordinatetransformationsfromthelocalto theglobalreferenceframesare

required. Now,thegeneralizedcoordinatevector

η ^g

^{containsallthe}coordinates ofvariouselementsandtheglobal

M ^g

^,

K ^g _M

^and

C ^g _M

^{areformedbyaddingeach}

termofthematricesofallniteelementsin thecorrectplace.

3.2.2 Flexible rotor model

Inthiswork,thecontinuouselasticrotorandthelaminationsaremodeledusing

thebeamnite elements,which arebasedonthe Timoshenkobeams,and the

rigid disks. Forboth elements, the analysis focuses onthe lateral vibrations,

andthereforetheaxialandtorsionaldegreesoffreedomarenotconsidered.

The rigid disks are prismatic (properties of all cross sections are equal),

homogeneous,straight,anduntwisted(theprincipalaxesofelasticityareequally

directed in space). Each disk assumesa lumped mass matrix, comprisesfour

DOF,andisanalogoustothe4-DOFrotormodeldescribedinsubsection3.1.2.

However, in order to unify the description of all the elements, the vector of

the nodal displacements, the generalized coordinates of the only node of the

element,is

η = [x, y, β x , β y ] ^T

^andthe

M

^,

G

^,

f

^aremodiedappropriatelyas

Thebeamelementsareprismaticandhomogeneousbeams,aspresentedin

Fig.3.4,thataccountforthesheardeformationinthebeamelementtransverse

(radial)direction.Specically,intheapproachusedthesheareectisincluded

in theshape functions (Chenand Gunter, 2005). Eachbeamelementhastwo

nodeslocatedattheendsofthebeam,andeachnodehasfourDOF.Thevector

ofthegeneralizedcoordinatesoftheelement,is

η = [x 1 , y 1 , β x,1 , β y,1 , x 2 , y 2 , β x,2 , β y,2 ] ^T .

^(3.21)

It describes the time-dependent displacements of the end-points of the nite

shaftelement. Withthisdenition,thegeneralizeddisplacementofaninternal

point of the element placed on the z axis (center line) can be expressed by

interpolation,usingEq.(3.13)as



where theshapefunction matrix

Φ

^{consists of}translationaland rotational in-terpolation shape functions of the Timoshenko beam. Each column of the

Φ

x z

y b y , 1

x 1

b x , 1

y 1 b y , 2

x 2

b x , 2

y 2

b y x b x y

Figure3.4: Beamelementandcoordinates

is an individual shape function vector, which represents the mode associated

with a unit displacement of oneof the coordinates of oneof the nodes.

Uti-lizingthederivationsfrom theliterature,theshapefunctions andenergiescan

be computed, according to Chen and Gunter (2005), for each nite element.

Then,asstatedearlier,theequationsofmotionfortheparticularelementscan

becomputed,and afterthat, the equationsofmotionfor thecompletesystem

areassembled.

Inthecompletesystem,wemayassumeonlythedisplacementatthenodes.

Thecomplete systemisexpressed in theglobal coordinate system;nowanew

displacement vector has all the node displacements of all the elements as its

components,namelyfor

P

^{nodes,eachwiththelocalj-th}displacementvectors

η _j

^,theglobaldisplacementvectoris

η ^g =

η ₁ ... η _P T

. Sucha

straightfor-wardtransformationispossiblebecausethedisplacementsinthecorresponding

axesin thelocalandglobalsystemshavethesamemagnitudesanddirections.

Thecompleterotorsystemisdescribedby

M ^g η ¨ ^g + (D ^g _M + ΩG ^g ) ˙ η ^g + K ^g η ^g = f ^g ,

^(3.23)

wheretheglobalmass,damping,gyroscopicandstinessmatricesareobtained

asexplainedearlier. Equation(3.23)featuresmanyDOF,itiscoupledandthus

notpracticalforcontrolengineeringpurposes.

3.2.3 Reduction of the number of degrees of freedom

In practicalcontrol problems, only a limitednumber of degrees of freedom is

usuallyrequired. ThenumberofDOFinthederivedmodelcanbereducedby

modaldecompositionandthentruncationofhigh-frequencymodes. Inthecase

of thesymmetricpositive denitematrices, themethod leadsto anuncoupled

systemin thefamiliar form of Eq. (3.13), where themeasurementvector

rep-resents the physical coordinates. In other words, the center-line of the rotor,

determinedbytheaforementionedmeasurementvector,isalinearcombination

themethodisasfollows.

First,the

M

^(where

M

^{equalstothenumberofDOF)}independentnormal mode shape functions (eigenvectors)and corresponding criticalspeeds

(eigen-values)areobtainedbysolvingtheM-DOFeigenvalueproblem. Weassumethe

free-freevibrationsoftheundampedsystem(Eq.3.23),which isconsideredfor

Ω = 0

^{. Thecritical speedsof thesystem, at which the}self-excited vibrations occur,can befoundbysolvingthefollowinggeneralizedeigenvalueproblem

K ^g − ω ² _k M ^g

φ ^m _k = 0,

^(3.24)

where

ω k

^,

φ ^m _k

^are

M

^{critical speedsand} correspondingmodeshapefunctions, respectively. Thesystem'sself-excitedvibrationsareequaltothesuperposition

oftheglobaldisplacementsas

η ^g _k (t) = φ ^m _k · cos(ω k t + ς k ),

^(3.25)

where

ς k

^{is the arbitrary phaseangle. In physical systems, the matrix}

K ^g

^is

symmetric, real, and positive semidenite (all of whose eigenvalues are

non-negative). Thematrix

M ^g

^{issymmetric,real,andpositivedenite. Therefore,}

critical speeds and mode shape functions are real. Furthermore, mode shape

functions,which correspondtodierentcriticalspeedsareorthogonal. In

con-trolengineering,itisacommonpracticetoobtainanorthonormalsetof

eigen-vectors(orthogonalandwithunitlength)forthemass-weightedstinessmatrix

(massmatrix becomes identity matrix). This mayoer certain advantages in

controlengineering, albeitthephysicalmeaningofparametersandvariablesis

lost. Hereweuseadierentprocedure.

Second,inordertopreservethephysicalmeaningofthesystem'sparameters

and variables we follow Lantto (1997), and we dene the matrixof the mode

shape functions asconsisting ofthe rigidbody mode and exiblemodeshape

functionsas

Φ ^m =

Φ ^m _rigid , Φ ^m _flex

. Theorderofthesystemisreducedbyusing

the

Φ ^m _flex

^{thatconsistsofthemodeshapefunctions}

φ ^m _k

^{,arrangedinthematrix}

accordingto theascending critical speeds

ω k ,

^{and thentruncated. In}

Φ ^m _flex = φ ^m ₅ ... φ ^m _4+M 2

,

^only

M 2

low-frequencyexiblemodeshapefunctionsare retained, from theoriginal

M

^{ones (also}zero-frequencyrigid modes, obtained from the eigenvalue problem, are dropped). We dene new orthogonal

zero-frequencyrigidbodymodeshapefunctions

Φ ^m _rigid =

φ ^m ₁ ... φ ^m ₄

y axis),tiltingmotion(aboutx andy axis),locationofk-thnode,andlocation

ofthecenterofgravityoftherotor. Inaddition,toaddsomephysicalmeaning

to thedisplacements thatresultfrom theexible modes, the

Φ ^m _flex

^{isscaled in}

such awaythat themaximal displacement(resulting from thesamemode)in

Now, it is possible to relate the physical displacement vector

η ^g

^{and the}

modalcoordinates

η ^m

^as

η ^g = Φ ^m η ^m .

^(3.27)

Equation(3.27)istheapproximationaftertruncationofhigh-frequencymodes,

henceitneglectshigh-frequencydynamics. Thesetofreducedanddecomposed

(coupledonlyby

G ^m

^{matrix)equationsofmotionforthestudiedsysteminthe}

modalcoordinatesandinthegeneralizedform are

M ^m η ¨ ^m + (D ^m _M + ΩG ^m ) ˙ η ^m + K ^m η ^m = f ^m ,

^(3.28)

wherethemodalmatricesandmodalvectorofforcesare

M ^m = ( Φ ^m ) ^T M ^g Φ ^m , D ^m _M = ( Φ ^m ) ^T D ^g _M Φ ^m ,

^(3.29)

G ^m = (Φ ^m ) ^T G ^g Φ ^m , K ^m = (Φ ^m ) ^T K ^g Φ ^m ,

^(3.30)

f ^m = (Φ ^m ) ^T f ^g .

^(3.31)

InEq.(3.28),thediagonalmatrices

M ^m

^,

K ^m

^andtheskew-symmetric

G ^m

^are

obtainedfrom theFEM model. However,thedampingmatrix

D ^m _M ,

^{which for}

thesystemdecouplingisassumedto bestiness-proportional,isobtainedfrom

themeasureddamping.

The damping factors and critical speeds of the few lowest bending modes

canbemeasuredusingmodalanalysis. Thedampingismeasuredinrelationto

thecriticaldamping(

ζ k = 1

^{). Therefore,thedampingelement}

d ^m _M,k

^inthek-th

DOFgoverningequationofform

m ^m _k η ¨ ^m _k + d ^m _M,k η ˙ _k ^m + k _M,k ^m η _k ^m = 0,

^(3.32)

whichisequivalentto

¨

η _k ^m + 2ζ k ω k η ˙ ^m _k + ω ² _k η _k ^m = 0,

^(3.33)

is determined as

d ^m _M,k = 2ζ k k ^m _M,k /ω k

^{. The damping factors for the}

higher-frequencymodes,whichcannotbemeasured,hastobeestimated. Usually,for

high-frequencymodes,

ζ

^{increases(uptothecriticalvalue)whenthefrequency}

isincreased. However,theestimationof

ζ

^isuncertain.

Equation (3.28) consists ofthe rigid body modesand the selectednumber

of distinctslightly damped oscillators, which representthe exible modes. As

previouslystated,

k

^{-thoscillatordescribesthefreevibrationofthemodeshape}

function

φ ^m _k

^{. The}correspondingoscillatorhascomplexpoles

s k = − σ k ± jω d,k =

− ζ k ω k ± jω k

p 1 − ζ _k ²

^{and,asitrepresentsthe}under-dampedcase,itsvibration decaysas

η k ∼ e ^{−σ k t} cos(ω d,k t + ς k )

^.

Usually, four to tenmodeswith thelowest criticalspeedsaresucient for

accuraterepresentationoftheexiblerotor. Thereducedmodeshapefunctions

matrixcanbecorrectedbymodalcorrection; alsotheFEM modelcanbe

iter-ativelycorrecteduntilthedierencesbetweentheexperimentalmodalanalysis

andthemodelaresmallenough.

Thepresentedmodelsassumethattheaxialmotionisindependentfromthe

0.0100 0.0465 0.1268 1.0650 0.2000 0.9400 0.5883 1.1850

Figure3.5: Sketchoftherotormodelwiththeimportantlocationsgivenin[m],

andtherigidbodymodeshapefunctions

Table3.1: CriticalspeedsobtainedfromtheFEMrotormodelareslightlyhigher

thanthemeasuredones. ThemeasurementsweredoneusingBrüer's&Kjeær's

modeanalyzer.

Mode

k

^{FEM model, Modalanalysis, Dierence Dampingratio}

ζ k

frequency[Hz] frequency[Hz] [%]

1 260.3 259.8 0.2 0.004118

2 539.0 526.7 2.3 0.002263

3 951.8 948.2 0.4 0.004345

3.2.4 Rotor of the test rig

In particular, the FEM representation of the rotor of the test rig comprises

P = 32

^{nodes,locatedontherotor'saxisofsymmetry(zaxis);this}corresponds to128-DOF system. Therotorshapewiththe locationsofitsmostimportant

nodes (beginning of the shaft, location of axial disk, sensors white circles,

radial actuators black triangles, end of the shaft) and the location of the

centerofthemass,ispresentedinFig.3.5.

Thedampingfactorsandcriticalspeedsof thethreelowestbendingmodes

oftherotorweremeasuredusingBrüer's&Kjeær'smodeanalyzer. Theresults

forthedampingandcriticalspeedsarepresentedin Table3.1. Thedampingof

theexible modesfor the complete rotoris higher thanfor the rotorwithout

couplings,aspresentedinPublicationIV,therstthreefree-freefrequenciesof

the shaftalone are 343, 653, 1154Hz. As the dampingof the high-frequency

modes, which arenotmeasured,isunknown,it isestimatedasin aworst-case

scenario(

ζ k>3 ≈ 0.001

^).

The rigid body modes and the rstthree exible modeshapefunctions of

therotormodel arepresentedinFig.3.5andFig.3.6,respectively. Ingeneral,

ascanbe seenin Fig.3.6, thenodes ofthe mode shapefunctions are closeto

the location of bearings and sensors. In Fig.3.7, the complete rotor, that is,

theshaft,radialAMB laminationsand couplings,underthemodalanalysis,is

Mode 1, f =260.3 Hz Mode 2, f =539.0 Hz Mode 3, f =951.8 Hz

end−A end−B

Figure 3.6: Sketch of the rotor model and the mode shapesof therst three

low-frequencyexiblemodes,in(y,z)plane,aredepicted. Themodeshapesare

multipliedby

0.2

^{inordertotintothesamegureastherotor. Thelocations}

oftheactuatorsandsensorsareindicatedwithredandgreenstars,respectively.

The rotoris atrest, i.e., there are noexternal forces, under free-freesupport,

and

Ω = 0

^.

Figure 3.7: Complete rotor, under the modal analysis, consists of the shaft,

axial disk, propellercoupling,AMB laminations ttedwith aluminumsleeves,

andadditionalaluminumsleevesfortheposition measurements.

In document Design and implementation of FPGA-based LQ control of active magnetic bearings (sivua 59-67)