Model Order Reduction Of Rotordynamical Systems

Preprints of the 18th IFAC World CongressMilano (Italy) August 28 - September 2, 2011can be replaced by that of a spring-damper device. One ofthe main characteristics of the resulting matrices is thatthey have coupling terms between the different degrees offreedom involved 1 leading generally to non-symmetricalsystem matrices Do and Ko .Such a modeling approach of the bearing results in asystem where the rotor and the housing are coupled toeach other via the stiffness and the damping matrices ofthe substitute spring-damper model. This is reflected inthe layout of the overall system matrices K and D asshown below 2 for a simplified model kkt12 t11 kt21 kt22 koK ko 0 000 kokh11 kokh210 ) ) rotor Kt 0 bearing )kh12 housing Khkh22(4)This exemplary stiffness matrix corresponds to a systemwith only two degrees of freedom for the rotor and forthe housing, which are coupled to each other through onedegree of freedom.2.4 Overall model of the machinePractically, it is not possible to generate each of the threesub-models separately. The FEM software tool delivers thefollowing model of the complete machineMz̈(t) (D G(Ω)) ż(t) Kz(t) Ju(t)y(t) Lz(t).(5)which consists of few hundred thousands of second-orderdifferential equations. It incorporates all of the three submodels and their dynamic interactions, however, in a coupled manner that does not allow a straightforward decoupling. Hence, there is a need for an automatic algorithmthat first partitions the system (5) into the sub-modelslisted above, allowing, if necessary, to reduce them separately, and then re-couple the resulting models to conducta complete system simulation (see Sec. 4 and 5).3. KRYLOV-BASED MODEL ORDER REDUCTIONH(s) y(t) Cx(t),(6)of order n N, where E, A Rn n , B Rn m ,C Rp n are matrices with constant coefficients, u(t) Rm , y(t) Rp , and x(t) Rn , respectively, the input,output and state vectors of the system. The transfermatrix of Eq. (6),H(s) C(sE A) 1 B.admits the infinite Taylor series’ expansion about s0(7)1 It is assumed here that the bearing does not move in the directionof the rotor’s axis.2 The layout of the matrix D is similar. C((A s0 E) 1 E)i (A s0 E) 1 B(s s0 )ii 0 Xmis0 (s s0 )ii 0The aim of Krylov-based model reduction is to find areduced model of order q n, whose first few momentsmis0 around a certain expansion point s0 match those ofthe original one (Freund (2003)). This approach is alsoknown as moment matching.A numerically efficient possibility to calculate such areduced-order model is applying a projection to the original model,ErArBrz } {z } {z } {WT EV ẋr (t) WT AV xr (t) WT B u(t),(8)y(t) {z}CV xr (t),Crby means of the so-called projection matrices V and W.For the choice of these matrices, the block Krylov subspaceKS , defined in e.g. (Grimme (1997)) is employed,KS (A, B) span{B, AB, · · · , AQ 1 B},where A Rn n , B Rn m , and S Q · m. Specifically,if the projection matrices are chosen such that Kq1 (A s0 E) 1 E, (A s0 E) 1 B colspan(V), Kq2 (A s0 E) T ET , (A s0 E) T CT colspan(W), then the first qm1 qp2 moments of the original andreduced models match, and the procedure is called a twosided Krylov method. In the so-called one-sided method,only one Krylov subspace is used with the common choiceW V and only qm1 or qp2 moments match (Grimme(1997)). Note that for the matrices V and W to haveappropriate dimensions, the order of the reduced systemshould be a common multiple of the number of inputs andoutputs.For the numerical computation of the matrices V and Wthe Lanczos, Arnoldi and two-sided Arnoldi algorithmsand their numerous improvements and modified versionsare used. For more details see e. g. (Antoulas (2005)) andthe references therein.4. PARTITIONING THE OVERALL MACHINEMODELConsider the linear time-invariant multi-input and multioutput systemEẋ(t) Ax(t) Bu(t), XAs already discussed in Sec. 2, the overall machine modelconsists of the model of the rotor coupled to that of thehousing by means of the equivalent spring-damper modelof the fluid bearings. Hence, the aim now is to separatethe rotor from the housing while preserving the completedynamical characteristics of the overall model.As it is free from coupling terms, the overall mass matrixM can be easily divided into a rotor and a housing part,Mt and Mh , respectively. It is also straightforward toassign the non-zero elements of the overall gyroscopicmatrix G(Ω) to that of the rotor Gt (Ω) as no gyroscopiceffects are present within the housing. The input andoutput matrices J, L could also be easily divided into twoparts corresponding to the rotor and the housing.3583

Preprints of the 18th IFAC World CongressMilano (Italy) August 28 - September 2, 2011FoBearing Model (Eq. 11)Rotor Model (Eq. 9).yr,o , yr,our yr.yh,o , yh,o4.2 Automation of the partitioning stepuh yh(reduced) Housing Model(Eq. 13)Bearing Model (Eq. 11)FoFig. 3. Schematic of the loop calculation involving thethree main sub-modelsThe decoupling of the overall stiffness and damping matrices K and D requires special attention as they expressthe coupling of the rotor to the housing through the fluidbearings. In order to preserve the overall dynamics, thematrices of the equivalent spring-damper model Do andKo are first calculated. These matrices are obtained bysubtracting the rotor and housing part Kt , Dt and Kh , Dhfrom the coupling terms in the K and D matrices, basedon the layout of Eq. (4). The input and output matrices ofeach of the three sub-models are then chosen so that thevelocity and position of the interface nodes at the rotor andthe housing can be extracted, so that the force interactionswith the fluid bearing can be calculated. Such a settingallows a loop calculation as shown in Figure 3.4.1 Mathematical model of the decoupled systemsNow, assuming that the overall model of the machine isavailable in the form (5), the decoupled systems of therotor and the housing are then, utMt z̈t (Dt Gt (Ω))żt Kt zt (Jt Jt,o ) Fo (9)ytLt zt ,yt,oLt,oand uhMh z̈h Dh żh Kh zh (Jh Jh,o )Fo (10)yhLh zh ,yh,oLh,owhere Jt,o and Jh,o are the newly inserted input matricesto account for the fluid bearing forces acting on the rotorand housing, respectively. Similarly, matrices Lt,o and Lh,oare the newly added output matrices to allow the forcesof the rotor and that of the housing acting on the fluidbearing to be calculated usingDo (ẏt,o ẏh,o ) Ko (yt,o yh,o ) Fo .(11)It is important to point out that the difference betweenthe new added outputs is included in Eq. (11), becausethe vector of the generated forces Fo is caused by therelative displacement and velocity of the spring-dampersystem representing the bearing. Caused by the fact thatthe velocity is needed, it is necessary to get over to a firstorder model (see Sec. 5). Please note also that the numberof inputs and outputs added to the sub-models dependson the interface between the bearing and the rotor fromone side and the bearing and the housing from the otherside. Hence, the number of inputs and outputs can be quitelarge once these interfaces involve a large number of nodes.As the models being considered are large scale models withseveral thousands of equations, it is almost impossible topartition their corresponding system matrices manually.Hence, there is a need for a numerically efficient algorithmthat automatically extracts the matrices of the three submodels once those of the overall model are given. This is infact possible, once the information about the correlationbetween equation number, node number and degrees offreedom is available. Fortunately, to the author’s knowlege,all FEM software tools can deliver such information. It isgenerally in a text file where each line contains the orderedsequence of node number, degree of freedom and equationnumber. This information can be expressed as a matrixdenoted as N Nn 2 , where N(i, 1) j and N(i, 2) krepresent the i-th equation, node number j and k-th degreeof freedom.The automation of the decoupling of a large-scale rotordynamical system which is given in the form of Eq. (5)with N is summarized in the following algorithm:Algorithm 1 Partitioning AlgorithmInput: M, D, G, K, J, L and N.1: Identify all equations involving a node which is affectedby a gyroscopic moment (rotor nodes) and merge theirnumbers in the setR {i [1, n] a, b [1, n] :N(i, 1) N(a, 1) G(a, b) 6 0}.2:Define the set of equations involving the housingnodes asH [1, n]\R.Knowing R and H, the mass, stiffness, gyroscopyand damping could be separately stored for rotorMt , K t , D t , Gt and housing Mh , K h , D h .4: Identify all elements K(i, j) 6 0 where i is out of Rand j is out of H. Store the identified stiffness elementsin the new bearing stiffness matrix Ko . Merge theequation numbers in the sets3:Ro {i R j H : K(i, j) 6 0}Ho {i H j R : K(i, j) 6 0}.Calculate the rotor stiffness matrix Kt and that of thehousing Kh by subtracting the bearing stiffness Kofrom K t and K h , respectively.6: Generate the new inputs Jt,o , Jh,o and outputs Lt,o ,Lh,o out of the equation numbers of Ro and Ho .7: Repeat step 4, 5 and 6 for the damping matrix.Output: The matrices of the three sub-models (9), (10),and (11).5:5. ORDER REDUCTION AND RE-COUPLING THEOVERALL SYSTEMNow that the three sub-models are on hand in a secondorder form, the last step before model reduction would beto transform them into a classical first-order state-spacerepresentation of the form Eq. (6). The following step is3584

Preprints of the 18th IFAC World CongressMilano (Italy) August 28 - September 2, 2011the same for the housing and the rotor and for simplicityonly shown for the latter.EtAt (Ω)ẋtxtz } { z } { z } z } { { I 0żt0Izt 0 MtKt (Dt Gt (Ω))z̈tżt(Bt Bt,o ) utz } { z } {0 0ut, Jt Jt,oFo! !ytLt 0ztyt,o Lt,o 0.żt0 Lt,oẏt,o {z } {z} yt CtCt,o(12)Fig. 4. Simplified FE-model of the gas compressor driveconsisting of: rotor (blue), bearing housing (green),the model of the fluid bearing (red ) and springdamper devices (gray).At this point, it should be noted that the inverse ofthe matrix Mt for large-scale systems should be avoided,which is why a matrix Et has been included in the systemrepresentation above. Moreover, even though it is not themost common way, this choice of Lt,o is needed to obtainẏt,o as an output to calculate the fluid bearing force (seeEq. (11)).The model of the rotor is always a parametric dynamical system with the gyroscopic matrix depending on therotational speed, as shown in Sec. 2. Thus, a parameterpreserving model reduction of this sub-model is reasonable, however, unnecessary for the class of machines considered in this work, where the model of a rotor normally consists of only 50–100 nodes, i. e. involving between200 to 600 state equations. The housing, in contrast, isa parameter-free system with more than 100,000 statevariables. Hence, it is advisable and advantageous to onlyreduce the housing model with a classical MOR methodand re-couple it to that of the rotor instead of reducingthe overall model with a parameter-preserving approach.The reduction of the housing sub-model by a two-sidedKrylov-subspace method as shown in Sec.3, leads to thefollowing reduced-order modelredEredh ẋh redAredh xh BredhBredh,o uhFo red ỹhChỹh ỹh,o xredh .Credh,o ỹh,o(13)One advantage of the partitioning/re-coupling approachpresented in this work is that it connects the sub-modelsthrough their inputs and outputs, which are quantitiesthat can not be eliminated by the MOR step, unlikethe state variables. A block diagram illustrating the recoupling of the overall system using the reduced housing isshown in Figure 3. To avoid the numerical and implementation disadvantages occurring from loop calculations, thereduced system Eq. (13) is re-coupled to the unmodifiedmodels of the rotor and bearing. This is done by combining Eqs. (9), (11) and (13). Nevertheless, the reductionintroduces a small error in ỹh , that propagates throughthe state vector x̃t to the output ỹt .The final result of partitioning, reducing and re-couplingis expressed by Eq. (14). Note that the inversion of thereduced matrix Eredh and that of Mt is numerically possibleand not time consuming, as both quadratic matrices arenon-singular and only have a few hundred rows.6. TECHNICAL EXAMPLETo illustrate the approach of this paper, a simplified FEmodel of the gas compressor drive shown in Figure 1 isconsidered. This machine, with a power of 16 MW, wasbuilt for an offshore platform in Norway. Its rotor hasa length of nearly 6 m and the housing has a quadraticbase with a length of 4 m. The simplification of the modelinvolves just considering the rotor, the fluid bearing, andtheir housings while approximating the rest of the housingusing several spring-damper-devices as shown in Figure 4.By taking the constraints of several degrees of freedominto account, the FE-modeling leads to a second ordersystem of the form Eq. (5) with 4204 states. The consideredmachine has two fluid bearings which are represented byspring-damper-devices and each assumed to be connectedto the rotor at a given node.The resulting system is first partitioned with the proposedalgorithm 1 within 1 s, leading to the rotor system with364 states and the simplified housing system with 3840states. Both systems’ terminal matrices are augmented by4 new inputs and outputs which allows the interactionforce transmitted through the bearings to be calculated.The resulting multi-input and multi-output system ofthe housing is then transformed into a first-order systemwith 7680 states and reduced down to 540 states by theKrylov-subspace method introduced in Section 3. Theexpansion point was chosen as s0 56 rad/s, leading toan excellent approximation of the transfer function withinthe bandwidth of interest (see Figure 5).In order to analyze smooth running of the machine, it isnecessary to observe the vibration amplitude at severalnodes of the housing for different rotor excitations. As therotational speed of the rotor constitutes the main vibrationsource, the bandwidth of interest for these excitations iscommonly chosen as twice the nominal rotational speedof the rotor. In the amplitude plot of a given housing3585

Preprints of the 18th IFAC World CongressMilano (Italy) August 28 - September 2, 2011 x̃ tẋredh At (Ω) 0M 1t Bt,oEredh 10 K o DoBredh,o Ko Do Ct,o0 Ct,o00Ct,o0Ct,o Eredh 0 1 Mt Bt0 Eredh 1 ut ,Bredhỹtỹhuh 0M 1t Bt,o 1Ct 00 Ko Do Credh,o redredAredh Bh,o Ko Do Ch,o 0Credh x̃txredh x̃t xred h (14)has then been applied, when necessary, to each of thesesystems allowing a drastic time saving in the simulation ofthe re-coupled overall system. However, a number of openquestions has still to be solved to optimize the presentedapproach, including model reduction of systems with alarge number of inputs and outputs, the consideration ofthe nonlinear effects in the equivalent model of the fluidbearings describing the dependency of the oil film on therotational speed, and the analysis of the error dynamicswithin the re-coupled system caused by the reduction step.REFERENCESFig. 5. Comparison of the original model (8408 states) andthe reduced model (1268 states); Vibration amplitudeof a housing node depending on different excitationfrequencies of the rotor node at a rotational speedΩ 3000 rpmnode presented in Figure 5, the nominal rotational speedis chosen to be around 3000 rpm, corresponding to an analyzed frequency interval of between 0 rad/s and 630 rad/s.The relative error resulting from reducing the model of thehousing and re-coupling it with the rotor is smaller than0.01 % up to 740 rad/s.The calculation of the frequency response of the original system of the housing for a single rotational speedtakes around 121 s with optimized tools, while that ofits corresponding reduced system only takes 10 s. A simulation of the complete re-coupled system including reduction requires around 24 s. Hence, for a critical-speedmap, i. e. a simulation with 15 different rotational speeds,the complete simulation with the original model wouldneed about half an hour while with the reduced model itwould take only 174 s including the decoupling, reduction,and re-coupling steps. This corresponds to a relative timesaving of over 90 % for this relatively small problem. Whenconsidering the whole machine shown in Figure 1, firstsimulations have shown that model reduction slashes thecalculation time for a critical-speed map from about threeweeks down to just a few hours!7. CONCLUSIONIn this paper, the problem of order reduction of largerotordynamic systems occurring in industry has beenconsidered. Under the common assumption that only thecomplete model of the machine can be obtained fromthe FEM software tool, an algorithm for partitioning thelarge drive system into the rotor, housing, and bearingsub-models has been introduced. Model order reductionAntoulas, A.C. (2005). Approximation of Large-ScaleDynamical Systems. SIAM, Philadelphia.Athavale, M.M. and Hendricks, R. (1996). A small pertubation cfd method for calculation of seal rotordynamiccoefficents. International Journal of Rotating Machinery, 2(3), 167–177.Daniel, L., Siong, O., Lee, K., and White, J. (2004).A multiparameter moment matching model reductionapproach for generating geometrically parameterizedinterconnect performance models. IEEE Trans. onComputer-Aided Design of Integrated Circuits and Systems, 5, 678–693.Freund, R.W. (2003). Model reduction methods based onKrylov subspaces. Acta Numerica, 12, 267–319.Friswell, M.I., Penny, J.E.T., Garvey, S.D., and Lees, A.W.(2010). Dynamics of Rotating Machines. CambridgeUniversity Press.Glienicke, J., Han, D.C., and Leonhard, M. (1980). Practical determination and use of bearing dynamic coefficients. Tribology International, 297–309.Grimme, E.J. (1997). Krylov Projection Methods for ModelReduction. PhD thesis, Dep. of Electrical Eng., Uni.Illinois at Urbana Champaign.Huebner, K.H., Dewhirst, D.L., Smith, D.E., and Byrom,T.G. (2001). The Finite Element Method for Engineers.Wiley-Interscience.Panzer, H., Mohring, J., Eid, R., and Lohmann, B. (2010).Parametric model order reduction by matrix interpolation. at – Automatisierungstechnik, 58(8), 475–484.Peng, L., Liu, F., Pileggi, L.T., and Nassif, S.R. (2005).Modeling interconnect variability using efficient parametric model order reduction. In Proc. of the Design,Automation and Test In Europe Conference and Exhibition. Munich, Germany.Yamamoto, T. and Ishida, Y. (2001). Linear and Nonlinear Rotordynamics: A Modern Treatment with Applications. Wiley-Interscience.Zienkiewicz, O.C., Taylor, R.L., and Zuh, J.Z. (2005). TheFinite Element Method: Its Basics and Fundamentals.Butterworth-Heinemann.3586

Model Order Reduction (MOR) techniques o er a solution to this problem, by making it possible to generate accurate reduced-order models that are able to approximate the behavior of the original large-scale ones, while being muc