Generalized modal reduction method for the dynamic analysis of rotating mechanical systems

The paper proposes modal reduction method of the dynamic systems composed of linear nonconservative subsystems coupled by nonlinear discrete couplings. Classical approach to the modal reduction is based on the transformation of the generalized coordinates by the real modal submatrix of the linear conservative part of the whole system. In case of modal synthesis method, transformation matrices are the real modal submatrices of the conservative part of mutually isolated subsystems. Rotating mechanical systems contain gyroscopic effects and other influences of rotation and damping. The paper introduces a generalized modal reduction method based on the complex modal values of the whole system or the isolated subsystems. Their complex eigenvalues and eigenvectors are used for transformation of the generalized coordinates and reduction of the number of degrees of freedom. The presented method is focused on vibrating rotating systems with gyroscopic and dissipative effects and nonlinear internal couplings. c © 2020 University of West Bohemia. All rights reserved.


Introduction
The rotating mechanical systems (e.g. high-speed gearboxes, bladed disks, rotors, turbochargers) are composed of many flexible and rigid bodies (below subsystems) mutually joined by flexible nonlinear discrete couplings. The mathematical models of these subsystems are nonconservative with nonsymmetrical matrices and after discretization by the finite element method have large number of degrees of freedom (DOF number). The standard numerical methods of dynamic analyses of the rotating systems with nonlinear couplings are very difficult to apply. On this account, different methods of dynamic modal reduction were developed. Standard methods are based on nodal coordinates reduction and reduction of the natural modes respected in dynamic response. The first is based on subdivision of the generalized coordinates into master and slave DOFs [5,7] which are computed directly from the master coordinates neglecting inertia and damping forces. Both standard methods in different modifications called substructuring or reduction techniques are the most widely used for the linear structures. A general knowledge of the substructuring techniques applied in commercial software is presented, e.g., in [1]. A comparison of various reduction techniques for flexible multibody dynamics was presented in [15]. One of the most suitable and well established methods for DOF reduction of large multi-body systems is the modal synthesis method [8,18,21]. The classical approach of the modal synthesis method is based on the reduction of the natural modes of conservative models of subsystems with respect to the dynamic response. This method was used at the author's workplace for vibration analysis of the screw compressors [17] and for modelling of rotating bladed disks [13]. Parametric optimization of gear drives from the steady-state vibration poit of view excited by the kinematic transmission errors was shown in paper [22]. The modal synthesis method was used for modelling of gear drives vibration influenced by the time dependend gear mesh [4,19]. An influence of the various level of DOF number reduction in rotor dynamics with flexible disks was studied in [23]. The coordinates transformation of all or of the chosen subsystems was used by means of their modal submatrices composed of low-frequency real eigenvectors of mutually isolated, undamped and nonrotating subsystems. This approach was generalized [20] by quasistatic consideration of frequency higher eigenvectors of subsystems, whereas DOF number of the reduced model was conserved as in the classical approach.
Rotatig mechanical systems contain gyroscopic effects and additional influences of rotation and dissipation [2,6,14,16]. On this account, the eigenvalues and the right and the left eigenvectors of such rotating subsystems are complex. All modes corresponding to the eigenvalues from the frequency spectrum of excitation must be respected in the reduced model. An accuracy of the reduced linear dynamic models can be tested by the orthogonality matrix test. The cross-orthogonality matrix can be used as an indicator of the accuracy of the DOFs number reduction [1,7]. The concept can be applied to the comparison of pairs of compatible eigenvectors of the reduced and full model using the modal assurance criterion (MAC).
The main aim of this paper is to present the generalized modal reduction method with reduction DOF number of the whole system or individual subsystems for modelling of the rotating multi-body systems with strong gyroscopic effects, damping and friction in couplings. This method is tested on rotating blade triplet linked by the rigid blade shrouding with friction in contact surfaces. An influence of different levels of DOF reduction on the dynamic response accuracy excited by the harmonic forces is shown. A comparison with the classical approach to the modal reduction is further discussed.

Modal synthesis method
Let us consider the mechanical system (rotor, blade packet, rings) which can be decomposed into N linearized rotating or nonrotating subsystems. In general, equations of motion of axially symmetrical subsystem j rotating with constant angular velocity ω 0 about own spin axis can be expressed in the matrix form [3,6] Vector q j (t) of the generalized coordinates has dimension n j (DOF number containing the generalized coordinates in the inertia or rotating frame). Mass matrix M j is symmetrical. Damping B j and stiffness K j matrices, including the external couplings with the frame, can be symmetrical or nonsymmetrical, ω 0 G j is the skew-symmetric gyroscopic matrix. Skewsymmetric circulatory matrix ω 0 C j in the case of modelling in the inertia frame is generated by a rotating internal damping and in the rotating frame by an external izotropic medium, respectively. The term ω 2 0 K j,ω is explicitly added for the centrifugal stiffening. In the case of modelling in the rotating frame, matrix K j,ω is reduced owing to softening under rotation. Skew-symmetric matrices ω 0 G j , ω 0 C j and symmetric matrix ω 2 0 K j,ω vanish for nonrotating subsystems. Force vector f C j expresses effects of internal linear or nonlinear couplings of the subsystem j with surrouding subsystems and f E j (t) is the time-dependent vector in which all excitation forces are listed, including the unbalance force.
The first step of modelling consists in the first-order formulation of the equations of motion (1) in the state space defined by state vector in the form [2,6,7] The equations corresponding to mathematical model (1) are where Let all modal values λ (j) Modal properties of subsystem j are expressed by the complex diagonal spectral matrix Λ j ∈ C 2n j ,2n j and complex couple right and left modal matrices U j ∈ C 2n j ,2n j , W j ∈ C 2n j ,2n j . These matrices satisfy the biorthonormality conditions [2,10,14] where E j is the identity matrix of the 2n j -th order. We chose for each subsystem j two sets of 2m j (m j ≤ n j ) so called master right and left natural modes corresponding to m j pairs of complex conjugate eigenvalues (diagonal elements of Λ j ) λ corresponding to master spectral submatrix State vectors u j in model (3) are transformed by the master right modal submatrices m U j ∈ C 2n j ,2m j mutually uncoupled subsystems into the modal coordinates as We note that for m j < n j , the frequency-higher natural modes usually contribute less to the subsystem vibration and their contribution in a dynamic response of the whole coupled system can be neglected. After modal transformation (8) and premultiplying of equations (3) by the transposed left master modal submatrices m W T j , with regard to the biorthonormality conditions (5), equations (3) becomeẋ Taking into account structure of state vectors u j , the eigenvectors of the subsystems can be written in the form The modal submatrices defined in (6) can be written as where are the right and left master modal submatrices of uncoupled subsystems in the original configuration space of generalized coordinates q j . Equations (9) can be rewritten in the forṁ The global form isẋ where Matrices m Λ, m R and vector x in (14) can be rewritten in the form Spectral submatrix m Λ includes the chosen eigenvalues λ ν of the all subsystems with positive imaginary part and left master modal submatrix m R includes corresponding eigenvectors r (j) ν . The complex conjugate eigenvalues are arranged in matrices m Λ * , m R * and the complex conjugate modal coordinates are arranged in vector x * . We can use the MATLAB built in ode45 solver for integration of submodel where Equations (17) are solved under zero initial conditions x(0) = 0 using ode45 MATLAB built-in solver. It incorporates Runge-Kutta method which is of medium order. According to (8) and (11), the vector q of generalized coordinates can be expressed as where right master modal submatrix m Q corresponds to m Λ. Model (17) of the coupled system in space of modal coordinates x (j) ν uncoupled rotating subsystems has m = N j=1 m j DOF number and for N j=1 m j < N j=1 n j is reduced. According to (2), (8) and (11), the dynamic response of the arbitrary subsystem j in the original generalized coordinates and velocities is real Complex modal coordinates x (j) ν of subvector are separated from vector

Modal reduction of the global model
If the linear part of elastic and viscous forces in couplings between subsystems can be excluded from the nonlinear couplings, vector f C (q,q) in (14) can be written in the form Global stiffness and damping matrices K C and B C describe the linearized forces in couplings, vector f F (q,q) expresses the nonlinear friction forces and q = [q T 1 , . . . , q T N ] T is the global vector of the generalized coordinates. All equations of motion (1) can be expressed in the global form where M , B, G, K, K ω , C are block diagonal matrices formed from corresponding matrices M j , B j , G j , K j , K j,ω , C j of the uncoupled subsystems. The first-order formulation in the state space is then where Chosen natural modes of the global model with linearized couplings are expressed by one pair of right and left master global submatrices corresponding to spectral submatrix Using transformation and biorthonormality conditions equation (25) can be rewritten in the configuration space of master modal coordinates of the system coupled by the linearized couplings. Its form iṡ where according to (27), is the left master modal submatrix. A dynamic response in generalized coordinates, according to (24), (27) and (29), is real Modal values λ ν , right eigenvectors T and corresponding complex conjugate modal values satisfy the definition relations and biorthonormality conditions (30).
Similarly to Section 2, the following submodel can be integrated instead of equations (32). Vector x and left master modal submatrix m R correspond to spectral submatrix m Λ = diag [λ 1 , . . . , λ m ] including eigenvalues λ ν = α ν + iβ ν with the positive imaginary part. We note that matrices N , P of the rotating mechanical systems are functions of the spin speed ω 0 . The calculation of modal values λ ν , u ν , w ν can be repeated for each value of the speed.

Application -modelling of rotating blades
The presented modal reduction method was tested using three rotating blades ( Fig. 1) with friction effects in contact surfaces of the blade shrouding. The flexible blades are fixed to a disk rotating with constant angular velocity ω 0 . The blades are discretized by FEM using 1D Rayleigh rotating beam elements derived in [9,13]. The position on the blades j (j = 1, 2, 3) is described in rotating coordinate system x j , y j , z j by three displacements u i and three small Euler's angles ϕ The rigid blade shrouds are fixed in end nodal points C j (j = 1, 2, 3) of the corresponding blades.
Mathematical model of the one decoupled rotating blade j with shroud (subscript B) can be written in the matrix form [3]   Let us consider a static preloading of the external blade shrouds by normal pressure forces N 0 (see Fig. 1). Consequently, normal forces N A , N B , radial R A , R B and axial A A , A B components of the friction forces act on the internal contact surfaces of the blade shrouds. All the contact forces are concentrated into central points of the contact surfaces in planes ξ A η A and ξ B η B , respectively. All contact forces acting on the central blade shroud are shown in Fig. 1. The resulting normal forces, acting in axes ζ A , ζ B perpendicular to contact surfaces, can be written in the form where N st is magnitude of the normal force resulting from equilibrium conditions of three blades loaded by pressure forces N 0 and centrifugal forces in case of constant angular velocity ω 0 of the disk. Contact stiffness k C is linearized [12] on the basis of contact pressure σ = Nst A ef , where A ef is the effective area of the contact surfaces. Vectors n A and n B of dimension 6 transform the blade displacements in nodal points C j into the normal displacements at contact points A, B.
Components of the friction forces acting on central blade shroud (see Fig. 1) can be expressed in dependence on slip velocities of contact surfaces in the form where f is the friction coefficient between contact surfaces of the blade shrouds experimentally and numerically investigated in [11]. The normal forces and components of friction forces acting on the outside blade shrouds in points A, B have opposite directions. Absolute values of the slip velocities are Vectors r A , r B and a A , a B of dimension 6 transform vectorsq C j of blade velocities in nodal points C j to the translational velocities in ξ A , ξ B radial directions and in η A , η B axial directions, respectively. The mathematical model of three mutually coupled blades (see Fig. 1) in general coordinates q = [q T 1 , q T 2 , q T 3 ] T loaded by pressure forces N 0 , centrifugal forces and external excitation can be expressed, according to (37)-(41), in the global matrix form where M , B, G, K s , K ω , K d are the block diagonal matrices having the structure Force vector f 0 expresses the loading by pressure forces N 0 and ω 2 T is the vector of centrifugal forces. A dominant significance has resulting vector f C (q,q) of all the contact forces which includes moreover linearized normal forces (38) and also nonlinear friction forces (39), (40) in couplings. This vector, after a transformation of the mentioned forces N A , N B , R A , R B , A A , A B into corresponding nodal points C j , can be expressed in the form The nonzero coordinates of vector f C (q,q) correspond to vectors q C j positions in global vector q. In the case of continuous contact between all blade shrouds, according to (38) up to (43), the resulting vector of contact forces can be rewritten in the form in short Linearized contact stiffness matrix K C in the reduced form (only with nonzero elements at the marked positions corresponding to vectors q C j in global vector q) is Finally, mathematical model (42) can be written in the form Let us find the solution of (47) in the form q(t) = q 0 +q(t) where constant vector q 0 satisfies the equilibrium conditions of the non-vibrating system

Vectorq(t) expresses the blades vibration satisfying the nonlinear equations of motion
where global stiffness matrix includes elastic properties of the blades and couplings, bending stiffening produced by centrifugal forces acting on the blade and softening as a consequence of the modelling in the rotating frame. Nonlinear vector f F (q,q) of the friction forces is defined as the third member on the right side in (44), where components of friction forces are expressed in (38) to (40) by means of q = q 0 +q andq =q.
Using the method presented in Section 3, state vector u, matrices N , P and vector p in model (25) have this structure According to (34), the blades dynamic response in space of generalized coordinatesq is where q ν , q * ν are pairs of the complex conjugate right eigenvectors.

Numerical experiments
The method presented in Section 3 and applied to rotating blades in Section 4 was verified using three rotating simple steel blades with shroud ( Fig. 2) described by parameters [13] listed in Table 1.
The contact stiffness was assumed to be k C = 5 · 10 6 N/m and the friction coefficient is varied from 0 to 0.2. Each blade was devided into N = 10 rotating beam finite elements using N nodal points. The rigid shroud is fixed to the last blade nodal point N = C j .  For illustration, several low-frequency eigenvalues of linear model of the three blades with smooth contact surfaces between blade shrouds (in (49) is f = 0 ⇒ f F (q,q) = 0) rotating by angular velocity ω 0 = πn 0 30 rad/s (n 0 is rotational speed in rpm) are presented in Table 2 including the description of the corresponding mode shapes. Bending stiffening effects under rotation increase imaginary part of eigenvalues. The improvement of the dynamic response using generalized modal reduction method compared with the classical approach to the modal reduction was analysed using reference reduced model of the three blades (see Fig. 1) with smooth contact surfaces between blade shrouds. A classical approach is based on transformation of the generalized coordinates in model (49) in the formq where m V ∈ R n,m is the master modal submatrix of the undamped (B = 0) and nonrotating (ω 0 = 0) blade triple with smooth contact surfaces (f = 0, f F = 0). Equation (49), using transformation (53) and orthonormality conditions can be rewritten into the form The reduced model (55) rewritten in the first-order formulation is characterized by matrices The difference of the several low-frequencyr eigenvalues λ (m) ν of the homogenous part of equation (55) with regard to exact eigenvalues λ ν of the homogenous part of equation (49) is shown in Table 3 for n 0 = 3 600 rpm and different numbers of master mode shapes. The cumulative relative error for 10 frequency lower eigenvalues λ (m) ν with positive imaginary parts for different numbers m of master mode shapes of the reduced model (55) is shown in Table 4.  To illustrate the applicability of the new modal reduction method (MRM), we consider fictitious harmonic excitation acting in y axis direction only on the central blade shroud in nodal point C 2 . Let us consider N st = 10 N and harmonic force F 0 cos ωt. Basic angular frequency ω of harmonic excitation is ω = ω 0 p S , where in practice p S denotes the number of the stator blades uniformly distributed around the circumference of the rotating disk. Corresponding excitation vector f E (t) in (49) and (51) is The nonzero element corresponds to generalized coordinate v (2) i for i = N in direction of axis y. As a reference excitation, we consider amplitude F 0 = 10 N and angular frequency ω = 1 131.5 rad/s corresponding to excitation frequency 180 Hz (blade rotating speed n 0 = 3 600 rpm, p S = 3) close to the first triplet eigenfrequencies of the linear model (see Table 2). Each blade was discretized by N = 10 finite elements. The original nonreduced model of the simple testing blade triplet had 180 DOF. This DOF number we can arbitrarily increase.
As an illustration, the time behaviour of central blade displacement v (2) N (see Fig. 1) in the end blade nodal point C 2 in y-axis direction of the reference blade triplet model (for m = 30 and Coulomb friction coefficient f = 0.2) is shown in Fig. 3 (top). Differences where nonzero elements correspond to generalized coordinates v (j) i and w (j) i for i = N and j = 1, 2, 3. We consider the excitation generated by force amplitudes F y = F z = 10 N with frequency 900 Hz which is close to the fifth eigenfrequency (see Table 2). The considered phase shift ωΔt 1 = π between blades 1 and 3 corresponds to p B = 60 (δ = 6 • ), p S = 15 (ω = 15ω 0 ) and blade rotating speed n 0 = 3 600 rpm. This operating state is characterized by large bending vibrations of the outer blades in opposite phase in xz plane. As an illustration, the time behaviour of the end blade nodal points displacements w   The dynamic response in modal coordinates, using the classical modal reduction method, is investigated by integration of m second order nonlinear equations (55) instead of the simple m first order nonlinear equations (35). Figs. 6-9 show the comparison of the time behaviour of tangential displacements w (1) N and w (2) N in time interval t ∈ 0; 0.05 s computed for harmonic excitation (59). Two variants for m = 10 and m = 30, n 0 = 3 600 rpm and f = 0.2 by both modal reduction methods (MRM) with the full model of the blade triplet are presented. As an illustration, a comparison of the computational times using the classical t C and the generalized t G modal reduction method indicates ratio τ (m) = t C (m) t G (m) . This ratio for presented case is τ (30) = 1.86.

Conclusions
The main objective of the paper is to present the new modal reduction method based on the complex modal values of the linear nonconservative part of the models.With regard to high computational costs, an application of this method is suitable especially for dynamic analysis of the large rotating systems with nonlinear couplings. The method enables dynamic analyses of the damped rotating systems including all rotation effects and nonlinear contact forces in internal couplings between subsystems. The classical approach is characterized by transformation of the generalized coordinates using the real modal submatrix of the linear part of the undamped and nonrotating system. The new approach is based on the transformation by the complex modal submatrix of the nonconservative linear part of the rotating system including all rotating and dissipative effects. The dynamic response in master modal coordinates is investigated by integration of the first order nonlinear equations whose number corresponds to identical number of the second order nonlinear equations, using a classical approach. Consideration of the chosen master complex mode shapes improves approximation of the damped gyroscopic structures behaviour in comparison with classical modal reduction in the basis of the real mode shapes of the undamped and nonrotating structures. In addition, the computational time is shorter. This fact was verified by means of numerical experiments with the rotating blade triplet with smooth contact surfaces between blade shrouds for eigenvalues and with friction in contact surfaces for the harmonically excited dynamic response.