Model reduction of rotor-foundation systems using the approximate invariant manifold method

This work presents a model reduction method suited for performing nonlinear dynamic analysis of high-dimensional rotor-foundation systems modeled by the finite element method. The approach consists in combining the component mode synthesis (CMS) method with the approximate invariant manifold method (AIMM), and allows the obtention of forced responses through the integration of a single pair of ordinary differential equations. The proposed approach is tested using two examples: a simple and a complex rotor-foundation system. In both cases, the nonlinearity comes from the fluid-film bearings. The results show that the method can provide a significant reduction in numerical cost while still retaining good accuracy when compared to direct time integrations. By means of the proposed method, the nonlinear dynamic analysis of high-dimensional rotor-foundation system becomes a feasible option.


Introduction
The complexity of rotating machines makes their design and analysis a challenging subject. One of several complicating factors is the interaction between the rotor and its supporting structure. This interaction has been shown to affect the critical speeds of rotor systems [1,2], which need to be avoided for a safe operation of the machines. Additionally, foundation flexibility is very relevant to the design of fluid-film bearings, which are commonly used in many rotating machinery applications [3]. These types of bearings are known to have unstable regions, denoted as "oil-whirl" and "oil-whip" [4,5], and the inclusion of the foundation flexibility is also known to affect the onset speed of this instability [6]. Aside from affecting the onset speed, the nonlinear dynamic behavior of the rotor when considering the foundation may be greatly affected during the oil-whip [7,8]. These works highlight the importance of including the flexibility of the foundation in rotor dynamic analyses.
The foundation of large rotating machines, such as turbine-generator systems for example, are often very complex structures [9,10]. For this reason, a classical approach is to introduce the foundation in the receptance matrix of the rotor system, after the equations have been transformed into the frequency domain [11][12][13][14]. In this case, experimentally obtained parameters of the foundation, such as modal masses, damping and natural frequencies, can be used, making the dynamic analysis in the frequency domain very accurate. These parameters can be obtained from run-down analysis [15] or operation modal analysis (OMA) [16]. This method, however, is not suited for nonlinear dynamic analysis in rotor-foundation systems, mainly because of the evaluation of the nonlinear forces. Thus, in these cases, the procedure is to model the foundation by means of the finite element method (FEM), and study the complete rotor-foundation system [17][18][19], which demands considerable computational power, and a reduction method is often required.
One well-known method used to reduce highdimensional models obtained by means of the FEM is the component mode synthesis (CMS), also known as substructuring [20,21]. The basis of the CMS is to divide the system into subcomponents or substructures. After this, a model reduction technique is applied in each subcomponent. The global system is obtained by assembling the substructures, leading to a reduced number of degrees of freedom (DOFs). Many different methods have been proposed to perform this subcomponent reduction and the subsequent assembly, such as the Craig-Bampton (CB) [22], MacNeal [23] and Rubin [24] methods. Substructuring is still an active field of research, and more information on the many different approaches can be consulted in the book by Allen et al. [21], and the reviews by Craig [20], Seshu [25] and de Klerk et al. [26].
Another approach to obtain reduced models is using the so-called nonlinear normal modes (NNMs), also known as invariant manifold method (IMM) [27,28]. The idea of the IMM is to extend the invariance property of eigenspaces of linear systems in the nonlinear range. This property states that any motion that starts on the eigenspace, remains on it as time tends to infinity. When nonlinearities are considered, however, the linear eigenspaces tend to be no longer invariant. The idea put forward by Shaw and Pierre [27] relies on constructing curved surfaces, or manifolds, that are tangent to the linear eigenspace and allows one to regain the invariance property for nonlinear systems. These are the invariant manifolds (IMs), which Shaw and Pierre also labeled as the NNMs of the system. However, the IMs obtained by the IMM are in general not unique, as discussed by Haller and Ponsioen [29]. The literature on IMs and NNMs is vast, and one may find more information in the reviews by Touzé et al. [30], Renson et al. [31], Mazzilli et al. [32], Avramov et al. [33] and Albu-Schäffer and Santina [34].
There are several ways to construct IMs for nonlinear systems. A common approach is to enforce a functional relation between a set of coordinates, labeled as master coordinates, of the system to another set, the slave coordinates. This functional relation, when substituted in the equations of motion, leads to nonlinear partial differential equations (PDEs), which solutions give the IMs. Some other ways to obtain the IMs are through shooting methods, harmonic balance or numerical continuation (see [31] for a review on the computation of IMs). Several methods have been developed to solve these PDEs. Taylor expansions are the most common way [35][36][37]. In this case, the relationships between the master and slave coordinates are expanded into polynomial series, and the coefficients of the polynomials are obtained from the PDEs. Here, one can also make use of the parameterization method [38], which is a powerful technique to obtain IMs. The three main parameterization styles are the graph style, the normal form style and the mixed style. In the first case, one has the already mentioned functional relationship between master and slave coordinates. In the normal form style, one has an additional nonlinear transformation in the equations of motion of the master modes, which is done in order to simplify them by removing resonant monomials and nonessential terms. This makes the method applicable to cases in which the IMs fold [39]. The mixed style is simply a combination of the two previously mentioned styles (one is referred to Touzé et al. [30], for more details). Alternatively, one may use numerical methods to solve for the IMs, such as the Galerkin method [40][41][42], FEM [43] or the finite difference method (FDM) [44]. Since the manifolds obtained from such numerical procedures are an approximation of the true IM of the system, one might label this approach as approximate invariant manifold method (AIMM).
This paper introduces a model reduction procedure applied to high-dimensional rotor-foundation systems. The method combines the use of CMS and the AIMM to obtain fast and accurate solutions when considering nonlinear bearing models. The CMS method used was the CB [22], which divides the rotor and foundation into the so-called superelements. The AIMM, based on [41], is then applied in the system reduced by the CMS, leading to a single pair of equations for the master coordinates. The main contribution of the work is the introduction of a method that allows nonlinear dynamic analysis to be performed in high-dimensional rotor-foundation systems.
The remaining of this paper is divided as follows: Sect. 2 presents the model of the rotor-foundation adopted. The source of nonlinearities comes from the fluid-film bearings, which are introduced in Sect. 3. The first reduction by the CMS is given in Sect. 4, and the subsequent application of the AIMM in Sect. 5. After this, Sect. 6 presents the results and the paper is finished by some conclusions in Sect. 7.

Rotor-foundation model
The rotor-foundation system is modeled by means of the finite element method (FEM). Figure 1 depicts the system considered, wherex r andx f denote the nodal displacements and rotations of the discretized mesh of both rotor and foundation, respectively. The form of these nodal vectors depends on the type of finite element utilized. For example, for the 3D finite element, one has, where u and ψ denote the displacement and rotation, respectively; and i denote the ith node in the mesh. Throughout this work, the bars denote that the displacements are not measured from the equilibrium position, which means that they denote the absolute displacements of the system. The equation of motion of the system can be written in general form as [13,14,19], where the subscripts r and f denote the rotor and foundation terms, respectively, M is the mass matrix, C is the damping matrix, G is the gyroscopic matrix, K is the stiffness matrix, and Ω is the shaft speed. The matrices R and S are the linearized damping and stiffness matrices of the bearings, while f nl are the purely nonlinear components of the bearing force. Also, the term f g denotes the gravity force, which is a constant vector acting at the center of mass of the system; and f h is the excitation due to mass unbalance, and it is given as, where f hy and f hz are horizontal and vertical Boolean vectors that define which node in the mesh the force is acting on (generally they act on the rigid disks), and m un is the unbalance moment. It is convenient to rewrite Eq. (1) with respect to the equilibrium position. This can be done by a change of variables as, being , and x e = {x er x ef } T the equilibrium position, which has the equilibrium displacements of both rotor (x er ) and foundation (x er ). The equilibrium position is obtained by assumingẍ r =ẍ f =ẋ r =ẋ f = f h = 0 in Eq. (1), which gives, Equation (4) is a nonlinear algebraic equation, and its solution provides the equilibrium position x e of the system. It is noted that the solution depends on the speed of the shaft Ω. By using (3), Eq. (1) is rewritten as, where, Also, it is worth noting that the inputs for the nonlinear force f nl are the relative displacements and velocities of the rotor and foundation.

Hydrodynamic bearing model
The bearings considered in this work are hydrodynamic. These types of bearings have many advantages, such as high load capacity and low friction, and are largely used in rotordynamics applications [7]. The dynamics of the fluid-film in the hydrodynamic bearing is commonly modeled using the Reynolds equation, which is given as [45,46], being p the pressure distribution, h the oil-film thickness, R the journal radius, μ the fluid viscosity (assumed constant), and θ and x denote the angular and axial coordinates. The bearing considered in this work is cylindrical, and its geometry is depicted in Fig. 2.
The fluid-film thickness is given by, with, where c r is the radial clearance, ε = e/c r is the dimensionless eccentricity, andū r y ,ū r z ,ū f y andū f z are the horizontal and vertical displacements of the rotor and foundation at the bearings. To obtain the bearing force, the Reynolds equation is solved based on the short bearing theory, which neglects the pressure gradient in the circumferential direction (∂ p/∂θ = 0). Thus, the pressure from the oil film can be obtained as [45], where L b is the bearing length andα denotes the whirl speed. In order to obtain the radial ( f r ) and tangential ( f t ) components of the fluid-film force, the pressure field has to be integrated over the axial and circumferential direction, that is, The integration limits θ 1 and θ 2 define the type of pressure distribution assumed. Here, the widely used Gümbel or half-Sommerfeld condition is used, where θ 1 = 0 and θ 2 = π . Thus, the integration above can be performed analytically, leading to, To express the force in cartesian components, one can apply the transformations (see Fig. 2), with, The expressions above model the hydrodynamic force in cylindrical bearings under the short bearing assumption, and it is widely used in the literature due to its simplicity [7,47,48].

Linearized bearing force
The linearized force can be obtained from a Taylor expansion of Eq. (14) around the equilibrium of the journal inside the bearing, that is, being Δu i and Δu i the relative displacements and velocities (i = y, z), k i j and c i j the stiffness and damping coefficients (i, j = yy, zz, yz, zy), and f y0 and f z0 the static forces. In the present case, f y0 = 0 and f z0 = F 0 , where F 0 denotes the vertical force on the bearings due to gravity alone. Note that, due to the change of coordinates (3), the equilibrium position is at the origin of the coordinate system. By performing the differentiations above, one arrives at [45,49], with, The eccentricity ε denotes the position of the journal inside the bearing, and it is dependent on the speed Ω, which in turn makes the coefficients above also dependent on Ω. One can solve the following nonlinear equation to obtain the eccentricity [45,49] The zeros of Eq. (19) give the locus of the journal center as the speed increases. The global bearing matrices R and S that are used in Eq. (5) can be obtained in the following way. Firstly, one can separate the rotor and foundation parts as, Secondly, one writes the sub-matrices S rr , S rf , and so on, as, being n b the number of bearings, and B r and B f Boolean matrices that indicate which rotor and foundation DOFs the bearings are acting, respectively. For an example, consider that the bearings act in the first node of both rotor and foundation meshes, and they are modeled using 3D finite elements. Recall that for 3D elements the ith , and the bearings act in the y and z direction. Thus, the Boolean matrices will be, In addition, the nonlinear part of the bearing force, f nl , can be obtained in a similar way, as where one subtracts the full hydrodynamic forces f y and f z given by Eq. (14) to the linear part f l y and f l z given by Eq. (16) to obtain the purely nonlinear contribution.

Component mode synthesis
The first reduction applied in the system given by Eq. (5) is performed using component mode synthesis (CMS). Since rotor-foundation systems modeled by the FEM can reach very high numbers of DOFs [17][18][19], the purpose of the CMS is to reduce this number for the application of the method presented in Sect. 5. The CMS approach selected was the Craig-Bampton (CB) method [21,22]. This method is well established in the literature and it is commonly used in the obtention of reduced order models (ROMs) for large systems modeled by the FEM. Some examples of the application of the CB method can be found in [50][51][52].
The basis of the CB method is the separation of the DOFs into internal and boundary or interface nodes, creating the so-called superelements. The boundary nodes are kept into physical form, while the internal nodes are reduced using fixed-interface modes. This makes the method very applicable to localized nonlinearities, such as friction interfaces [53]. In the present case, the boundary nodes correspond to the DOFs where the bearings are located at both the shaft and the foundation. Thus, Eq. (5) is rewritten as, where i and b denote the internal and boundary DOFs. The relation between the boundary and the internal DOFs is obtained by performing a static reduction, which gives the following reduction basis, The basis ψ ψ ψ has a size N × n B , being N the total number of DOFs in Eq. (24) and n B the number of boundary DOFs. Equation (25) is obtained by performing a unit displacement in each of the boundary DOFs and fixing all remaining DOFs. In order to improve the dynamic analysis, the static reduction is augmented by the use of the vibrating modes of the structure with all boundary DOFs fixed. These modes are obtained from the solution of the following eigenvalue problem: where n I denotes the number of modes retained. The eigenvectors φ φ φ r and eigenvalues λ C M S r correspond to the vibrating modes of the system when the boundaries are considered fixed. The basis with the fixed vibrating modes will be given as, Thus, the complete CMS basis is obtained as, which has size N × (n B + n I ). Using the CMS basis, the displacements are written as, being q i the modal (generalized) coordinates of the fixed-interface modes. By substituting the expansion (29) into Eq. (24), and multiplying by T T , one obtains the reduced set of equations for the rotor and foundation, that is, where, Equation (30) is the equation of motion of the rotorfoundation system after the reduction into superelements. Thus, q i consists of the internal coordinates of both the rotor and foundation, and x b the boundary coordinates at the bearings, which is where the rotor is connected to the supporting structure. Figure 3 depicts the CMS reduction process in the rotor-foundation system.

Approximate invariant manifold method
The next step in the present approach is the application of the approximate invariant manifold method (AIMM). The basic idea of the AIMM is to estimate the IM of the dynamical system, allowing one to obtain the response of it by only solving a selected set of master coordinates. To this end, the equations reduced by the CMS are projected into the eigenspace, that is, the space spanned by the vibrating modes of the coupled structure.
Here, it is important to note that one could bypass the application of the CMS presented in Sect. 4 and apply the modal analysis directly into the full system given by Eq. (5). However, if the model has a very high number of DOFs (which is common in systems modeled by the FEM), the solution of the eigenvalue problem may be too computationally expensive. Hence, the use of CMS can be very useful in such cases [21]. Alternatively, the computational burden can be avoided by using the direct computation of IMs [37,39], which only requires the computation of the master modes instead of all the eigenvectors of the system.
Due to the damping of the bearings and gyroscopic effect of the shaft, a complex modal analysis is necessary to fully decouple the equations at linear order [54].
One first writes the equations in the state-space form as, where, Next, the state vector is expanded in terms of the eigenvectors of the matrix A as follows, where n are the modes retained, η η η is the modal matrix of size 2N × 2n that has the eigenvectors η η η i on its columns, and q(t) are the modal (generalized) coordinates. Here, one has the opportunity to perform a second reduction in the equations. In the CMS, the system is reduced from N to (n B + n I ). In some cases, specially in highly discretized 3D finite element models, the number of boundary DOFs n B may still be large. Thus, one can further reduce the system by choosing n in Eq. (34), such that 2n < (n B +n I ). The eigenvectors and adjoints are the solution of, being λ i the eigenvalues, which are generally complex conjugate pairs,η η η i the ith adjoint eigenvector, and H and * denote the hermitian (complex conjugate) transpose and complex conjugation, respectively. The adjoint eigenvectors are the complex conjugate of the left eigenvectors, and they are necessary here due to the non-symmetric nature of A, which in turn is due to the gyroscopic effect and anisotropy of the bearings. The relation between the eigenvectors and its adjoint is given as, being δ i j the Kronecker delta. Note that the previous relations assumes that the vectors have been normalized. By substituting the expansion (34) in the equation of motion (32), and multiplying by the adjoint matrix η η η H , one has, which are 2n, first-order, complex equations uncoupled at linear order. By separating the modal coordinates into its real and imaginary parts as q i = p 2i−1 + j p 2i , with j = √ −1, Eq. (37) can be written using real numbers. The result is, with σ i = Re{λ i } and ω i = Im{λ i }. Note that the procedure above assumed that the mode i is a vibrating mode, that is ω i = 0. In case, ω i = 0, one has a ith overdamped mode, and its equation of motion is simply, where the eigenvectorη η η i is now purely real as well. The solution manifold corresponding to Eq. (38) is 2D while the ones given by Eq. (39) is 1D [29].
To apply the AIMM, one chooses one mode, or a pair of coordinates, from the 2n used to express the displacements in Eq. (34). These are then used as a basis to describe the remaining 2n − 2 slave modes of the system. It is worth mentioning that the approach could be easily extended to multiple master modes, and there is no loss of generality here. Let u = p 2k−1 and v = p 2k denote the selected master coordinates. The relation between the remaining modes with the master coordinates can be expressed as, where P i and Q i can be seen as the coordinates that compose the ith manifold. Since the system is non- Fig. 4 In the direct numerical integration, the solution is obtained as a trajectory in phase space (a). In the AIMM (b), one obtains the trajectory of the master mode only, and uses the shape of the manifold to obtain the full trajectories autonomous (forced) due to the unbalance force, one needs to introduce a third generalized coordinate φ = Ωt, which is the phase of the external excitation. Since P i and Q i depend on three quantities, one can think of the manifold as a 2D surface with time-varying coordinates due to the external excitation [41]. The procedure followed in the AIMM is depicted in Fig. 4. In the traditional direct integration, the solution of Eq. (38) gives the trajectories in phase space (depicted in blue). In the AIMM, instead of integrating the full equations, one only solves for the master coordinates u and v, and uses the relations in Eq. (40) to obtain the trajectories of the full system. By differentiating Eq. (40) with respect to time one arrives at, where the functional dependence of the variables has been omitted for better clarity. Substituting Eq. (38) above, and consideringφ = Ω, one has, where . Note that the external vector G depends now only on u, v and φ due to Eq. (40). The above equations are partial differential equations (PDEs), and they can only be solved approximately by means of numerical methods. Their solution gives the coordinates of the manifolds P i and Q i , for the ith mode. In the AIMM, the coordinates are assumed as, being C lmr i and D lmr i the unknown coefficients, U l,m and S r known shape functions, and N u , N v and N φ the number of shape functions assumed. For the expansion in u and v, the shape functions were assumed standard Chebyshev polynomials, which are known to be very accurate in a wide range of applications [55,56]. In this case, the 2D shape functions are obtained through the tensor product of two 1D polynomials in the u and v directions. For the expansion in φ, a Fourier series expansion was performed, taking advantage of the periodicity of this coordinate [41].
The coefficients of Eqs. (43) and (44) are obtained from a weighted residual method, namely the Galerkin method. Firstly, the form assumed in Eqs. (43) and (44) is substituted in (42), leading to, The next step is to perform a Galerkin projection. Thus, one multiplies the residue above by the same shape functions and integrate the result, that is, for i = 1, 2, . . . , n − 1; p = 1, 2, . . . , N u ; where R 1i and R 2i are the residues of Eqs. (45) and (46), respectively, and [u 1 , u 2 ] and [v 1 , v 2 ] are the integration limits; they denote the region of validity, in the phase space, of the approximate manifolds. The integrations in Eq. (47) were performed using the roots of the Chebyshev polynomials in U p,q (u, v), also known as Gauss points [55]. The number of points used were N v + 1 and N u + 1. In addition, in the integration for the φ variable, the number of points must be at least 2N φ in order to minimize aliasing errors [56]. Equation (47) consists of integro-algebraic equations that need to be solved for the coefficients C lmr i and D lmr i . The total number of equations will be N u N v N φ (2n − 2). One may note that this number may be very high, and thus the method is useful mainly when the system can be described with a small number of vibrating modes n. It is also for this reason that a spectral method, instead of a finite difference or finite element approach, was chosen for the solution of (42). Spectral methods allows accurate solutions for boundary-value problems with less computational costs [55].
Since the idea of the AIMM is to find an approximation of the manifolds defined in (40), there may exist multiple solutions to the Galerkin equations (47). In order to obtain fast and accurate solutions, one is advised to use the solution of the underlying linear system as an initial condition for the nonlinear solver. This linear solution is obtained by neglecting the nonlinearities of the bearings in Eqs. (45)- (47). The integration limits [u 1 , u 2 ] and [v 1 , v 2 ] can also be estimated based on the solution of the linear system.
With the functions P i and Q i , for i = 1, 2, . . . , n, at hand, the equations of motion for the master coordinates can be solved, which are given by, The physical displacements and velocities are then obtained by, To summarize the present approach: the full equation, Eq. (1), obtained after finite elements discretization, is firstly reduced using CMS, Eq. (29). This reduced set is then written in terms of the normal modes of the structure (and possibly further reduced), Eq. (34). Lastly, the equation of motion for the modal coordinates is reduced to a single master mode, Eq. (48). Therefore, the proposed method reduces a N degree of freedom system, with N 1, to a single pair of equations. Provided the solutions of the manifolds are precise, Eq. (47), the AIMM can provide very accurate and fast solutions for high-dimensional dynamical systems.

Results and discussion
This section presents applications of the method proposed in the previous sections. The method will be studied in two systems: a simple rotor on a springmass foundation, and a realistic model of a turbomachine on a plate-like elastic foundation. In all sub-sequent simulations, the results were obtained using MATLAB TM . Numerical integrations were performed using the ode15s integrator, which is ideal for numerically stiff systems [57], with standard options and zero initial conditions. The Galerkin equations (47) were solved by means of the fsolve function, with the initial conditions being the solution of the underlying linear system.

Simple rotor-foundation system
This first example is dedicated to show how the proposed method works, and hence a simple system was chosen for the study. The rotor system is depicted in Fig. 5, and it consists of a shaft with a disk positioned at its midspan and two identical bearings at its extremities. In contact with the bearings is also the foundation, which is considered to be two spring-mass-damper systems, as depicted in the figure. The supporting stiffness and damping are considered isotropic. The shaft is modeled by Timoshenko beam elements, in both the y and z directions. Only bending motion is considered, thus torsion and axial displacements are neglected. The mesh has a total of 24 elements and 25 nodes, each with 4 DOFs (two displacements and two rotations in y and z). All the relevant data of the model is listed in Table 1.
The total number of DOFs in the system is N = 104. The first reduction consists in applying the CMS in the rotor. The boundary DOFs are the displacements and rotations of nodes 1 and 25 (8 in total) together with the foundation displacements (4 in total), while the remaining DOFs are labeled as internal. To construct the CMS basis given by Eq. (28), the first six fixedinterface vibrating modes were considered in the matrix φ φ φ. This leads to a reduction from N = 104 to n B +n I = 18, being n B = 12 and n I = 6. Next, the AIMM is applied to the equations reduced by means of the CMS. The number of modes retained in the modal expansion in Eq. (34) was 2n = 20. From these retained modes, 4 are highly overdamped modes, while 16 are vibrating modes coming in complex conjugate pairs. The inclusion of the overdamped modes makes the modal equations, Eq. (37), stiff, but they are necessary for the nonlinear analysis. Figure 6 shows the eigenvalues of the 6 vibrating modes considered, where FW and BW denote forward (shaft whirl in the same direction as rotation) and backward (shaft whirl in the opposite direction as rotation) modes, respectively. A distinct characteristic of rotating shafts is that the natural frequencies change with the increase in speed, due to the gyroscopic effect and the hydrodynamic bearings. The line ω = Ω in Fig. 6a indicates where the speed equals the natural frequencies and are denoted the critical speeds of the system [45,49]. The first backward and forward critical speeds are found to be ω B c = 1092.2 rpm and ω F c = 1102.1 rpm (Fig. 6a). A crucial step in the application of the AIMM is the selection of a master mode to enslave the remaining modes of the system. The procedure to follow here is to choose the slowest mode, that is, the mode with the smallest absolute value of the real part of the eigenvalue. This is frequently called slow-manifold reduction [29], and it is a common choice in model reduction. From Fig. 6b, it is clear that the 2 FW mode has the Fig. 6 Eigenvalues of the simple rotor-foundation system: imaginary (a) and real part (b) Fig. 7 Comparison between the reference solution with the AIMM considering different number of shape functions for Ω = ω F c and m un = 0.05 kg · mm: radial displacement (a) and orbit (b) of the rotor at the bearing smallest real part up to around 2000 rpm, where the 2 BW mode actually becomes positive, making the system unstable. This phenomenon is the so-called "oilwhip" and is very well documented in the literature [4,5]. Since this instability is commonly avoided in the operation of actual rotating machines, the analysis was restricted to speeds below 2000 rpm and the master mode is considered to be the 2 FW. Another approach would be to choose two modes to serve as master modes, in this case 2 FW and 2 BW. This is known as multi-mode invariant manifolds [58]. This is not performed here, however, as the present approach only allows one master mode to enslave the remaining ones.
With the master mode chosen, one needs to set the number of shape functions used in the expansions (43)- (44). In order to define these numbers, convergence tests are necessary. Figure 7 presents one of the tests performed, where the displacements at the bearing are shown. The reference solution corresponds to the integration of Eq. (30). This result is obtained for the rotor at the critical speed, Ω = ω F c , and with an unbalance moment of m un = 0.05 kg · mm. As one can note from the figure, a convergence is reached for N u = N v = 3, which correspond to quadratic base, and N φ = 5, which correspond to four harmonics plus the constant term in the Fourier series. Therefore, these numbers are used in the following analysis. Figures 8, 9 and 10 show the radial displacements and orbits of the rotor at the bearings and disk positions and the displacements of the foundation, respectively, for different unbalance moments and for a speed close to the critical speed Ω = 0.9ω F c . Due to the symmetry of the system (see Fig. 5), only the displacement of one bearing and foundation spring-mass system is shown. The linear solution is also shown to assess the degree of nonlinearity. For m un = 0.05 kg · mm, one sees in Figs. 8a, 9a and 10a that the difference between the nonlinear response and linear is small, indicating a weak nonlinearity. As the unbalance is increased, this difference grows, and for m un = 0.5 kg · mm, the linear response gives discrepant results. One also notes that the nonlinearity is stronger in the bearings and the major effect of the nonlinear fluid-film force is to distort the orbits, which are ellipses in the linear case. By comparing the reference solutions with the AIMM, one notes good agreement, in both rotor and foundation responses, even at highly nonlinear cases. To obtain the reference solutions, 2(n B + n I ) = 36 equations are numerically integrated (as the system needs to be in state-space form to apply the numerical integrator), while in the AIMM only 2 equations are solved. Therefore, the method provides a great reduction while giving good accuracy in the responses. Figures 11 and 12 show the approximate manifold for the first (first overdamped) and fifth (first vibrating) modes in the case of strong nonlinearity. These manifolds are obtained from the solution of the Galerkin equations (47), and need to be obtained prior to the solution of the master equations of motion (48). Indeed, it is the correct obtention of these manifolds that allows good accuracy in the responses presented previously. As one notes from these figures, the manifolds are curved surfaces that move about with the phase of the unbalance φ. Also, this motion is not simple, as shown mainly in Fig. 11, where beside the translation of the surface, one notes additional "wobbles" (i.e., the shape of the surface changes with time), as described in [59]. In spite of this, it is clear that the translating motion of the surfaces has a more prominent effect than the change in their shape, which indicates that the expansion in φ bears a higher degree of importance in the obtention of the manifolds.
It is also worth mentioning that the AIMM only gives the steady-state solutions of the system. This fact makes the AIMM ideal for the obtention of steady-state solutions of rotors with low damping, where the transients take in general very long to die out, requiring a long simulation time. Figure 13 illustrate this by showing the full simulation time of the displacement at the bearing for Ω = 2000 rpm and m un = 0.1 kg · mm. It is clear that only a single cycle is enough for the obtention of steady-state solution with the AIMM, while in the full system, the whole simulation time is required, which in this case was 400 cycles.
The computation time required to obtain the results in Figs. 8, 9 and 10 is listed in Table 2. The application of the AIMM is in general more costly than the direct integration when one takes into account the solution for the manifolds. The numerical integration of the master equations, however, is very cheap, since they are only two equations, and it takes only seconds to complete.

Realistic rotor-foundation system
In this example, in order to evaluate the applicability of the proposed method, a more complex rotor system is studied. The mesh of the rotor is shown in Fig. 14a, and it consist of a multi-stepped shaft with three disks and two bearings. The rotor has a total length of L = 2.5 m, weights W r = 17.94 kN, and it is discretized into 26 elements and 27 nodes, where two 1D beams are used in both orthogonal directions (torsional and axial motions are neglected). Thus, each node has 4 DOFs, being two translations and two rotations. The mesh data can be consulted in [60] or [61]. The rotor is connected to the foundation through the bearings, positioned at nodes 6 (bearing 1) and 23 (bearing 2). The bearings are cylindrical with diameters d b1 = 160 mm and d b2 = 180 mm, lengths L b1 = 88 mm and L b2 = 98 mm, radial clearances c r 1 = 0.12 mm and c r 2 = 0.135 mm and fluid viscosity μ = 0.027 Pa · s.
The foundation consist of a steel plate of size 120×1200×3000 mm 3 , weighting W f = 33.27 kN. It is modeled using 3D linear hexahedral elements, with 8 nodes per element and 3 DOFs per node (three translations in x, y and z) [62]. The mesh consists of 208 elements and 350 nodes, and it is shown in Fig. 14b. The rotor is connected to two nodes of the plate at the top surface, as shown in the figure. Also, the plate is held in place by isotropic supports with stiffness K f = 10 10 N/m and damping C f = 5.8 × 10 5 Ns/m at its four lower surface vertices. The material properties of the rotor and foundation are the same as listed in Table 1.  The total number of DOFs in the system is N = 1158. Prior to the application of the AIMM, the system is reduced by means of the CMS. The boundary DOFs in the rotor are the nodes with bearings (8 DOFs), while in the foundation they are the DOFs with elastic supports (12 DOFs) and the connection with the rotor (6 DOFs). The number of fixed-interface modes retained for the rotor and foundation are 6 and 4, respec-tively. These numbers were obtained through nonlinear dynamic analysis, and provided satisfactory results. With the CMS, the system is reduced to n B + n I = 36 DOFs, with n B = 26 and n I = 10. Figure 15 shows a comparison between the first 7 eigenvalues of the full system and the reduced one. One can note good accuracy with the reduced system up to mode 6. Figure 16 shows the first 6 vibrating modes of the rotor-foundation system at Ω = 4000 rpm, where the displacements were rescaled for better clarity. Also, only the displacements of the upper surface of the foundation is shown. From these figures, one notes that the first two modes, (a)-(b), correspond to the conical and cylindrical modes of the rotor, and the fourth and fifth, (d)-(e), the elastic modes. The third (c) and sixth (f) modes are foundation-dominant, and, in the latter mode, the amplitude of the foundation is much higher than that of the rotor, as its displacement is weakly seen.
The next step consist in applying the AIMM to the system reduced by the CMS. Firstly, one needs to choose a master mode to enslave the remaining modes.
The best strategy is to find the mode with the lowest absolute value in the real part in the eigenvalues [29]. By looking at Fig. 15b, it is clear that the best candidate is mode 6, which is the slowest mode up to around 5000 rpm, where the system becomes unstable due to oil-whip. Since this instability is not considered here, the master mode chosen is the sixth mode. Here, it is also worth mentioning that the rotor-foundation system has 8 overdamped modes. As discussed in the previous example, these modes must be included in the nonlinear dynamic analysis. Due to the complexity of the system, no further reduction, without losing accuracy, was possible in Eq. (34). Thus, the number of modes Fig. 16 First six vibrating modes of the rotor-foundation system at Ω = 4000 rpm (Note: the displacements of both rotor and foundation were rescaled by 1000× for the sake of clarity, and only the upper surface of the plate is shown) considered was 2n = 2(n B +n I ) = 72, being 8 of these overdamped and 32 vibrating modes. Additionally, the number of shape functions used in expansions (43)- (44) where N u = N v = 3 and N φ = 5, which were obtained from a convergence study similar to the previous example (Fig. 7).
From the Campbell diagram, Fig. 15, one notes that a safe operating range for the rotor system, that is, with no critical speeds, is between 3000 and 4000 rpm. Thus, this was the range studied. In addition, the unbalance of the real turbomachine in which the model is based upon is estimated as m un = 0.0213 kg · m, and it is placed at disk 2 (node 15 in Fig. 1a). Figures 17, 18 and 19 show the radial displacements and orbits at several points of the rotor-foundation system for Ω = 3000 rpm and Ω = 4000 rpm: of the rotor at bearing 1 (node 6 in Fig. 1a), disk 2 (node 15 in Fig. 1a), and of the foundation at the connection with bearing 1 (Fig. 1b). Similarly to the previous example, the linear solution is also shown to illustrate the degree of nonlinearity. For Ω = 3000 rpm, the nonlinear effect is weak, and the linear response is very similar to the nonlinear one. At this speed, the highest excited mode corresponds to a foundation-dominant mode, namely mode 3 (see Fig. 16), and this is seen in the response by the higher vertical than horizontal displacements (due to the foundation bending). This highlights the importance of considering the flexibility of the foundation in the response. When Ω = 4000 rpm, the nonlinearity is apparent, specially at bearing 1, Fig. 17c-d. The main mode excited here is mode 6. By comparing the results of the reference solution with the AIMM ones, one notes good agreement. The AIMM provides accurate solutions for the high-dimensional system by the integration of 2 equations, leading to a great reduction in numerical cost, as the total number of DOFs of the system is N = 1158. Figures 20 and 21 show the approximate manifolds for the first overdamped and first vibrating modes of the system for different unbalance phases φ. These manifolds are obtained from the solution of the Galerkin equations (47). There are a total of 70 manifolds, which correspond to the slave modes of the system. Similar to the previous example, one notes that the manifolds are curved surfaces in motion. This motion, and the shape of the surfaces, is what allows a good agreement between the solution of the pair of equations (48) and the full high-dimensional system.  number of slave modes, the solution of the manifolds is computationally expensive and often takes longer than direct numerical integration. One can reduce the number of slave modes, sacrificing some accuracy, to decrease computational cost. It is important to note that integrating the master equations is fast as they consist of only a pair of equations.

Discussion
It is worth mentioning here some numerical aspects of the AIMM. The main disadvantage of the method lies in the solution of the manifold equations, Eq. (42), which are highly nonlinear PDEs. Different methods were employed in the literature to solve these equa-  Results obtained with a laptop with an Intel(R) Core(TM) i7-7500U CPU @ 2.90 GHz

Fig. 19
Radial displacements and orbits of the foundation at connection with bearing 1 for: Ω = 3000 rpm (a-b) and Ω = 4000 rpm (c-d) tions [31]. Here, the Galerkin method is used, which presented accurate solutions. However, this approach still requires a large computational cost, as the number of equations to be solved are (2n − 2)N u N v N φ , making the method applicable mainly to systems that can be reduced with a small number of modes. In the first example, since the system was simple, the solution of the manifold equations was very fast as shown in Table 2. In the second example, with the more complex rotor system, the solution took a bit longer, since the number of equations was larger, as Table 3 shows. Moreover, when the initial conditions are obtained from Approximate invariant manifold P 9 (u, v, φ) (first vibrating mode) of the realistic rotor-foundation for Ω = 4000 rpm and varying phase: φ = 0 (a), φ = π/4 (b) and φ = π/2 (c) the underlying linear system, the convergence is very fast, taking 4 to 5 iterations using the fsolve function of MATLAB TM . One should bear in mind that the numerical cost associated with the solution of the manifold equations is an offline cost, which means that they do not affect the actual solution of the system. Similar numerical costs exist for reduction methods based on the proper orthogonal decomposition (POD), for example, which is a well-known method to obtain reduced order models [63,64]. Given that the manifolds were solved, the reduction provided by the AIMM is simply astonishing: with just two equations the response of the whole system is accurately predicted. Not only this, but the solution gives the steady-state responses directly (see Fig. 13). Thus, since the cost of numerical integration is negligible, the only cost in the AIMM is the solutions of Eqs. (42).
The nonlinearity considered in this work was due to the fluid-film bearings. However, the method proposed is very general, and can be applied in a wide range of nonlinearities. Some possible applications are rotor-stator rubbing [65], geometric nonlinearities from large displacements [66,67], forces from seals [68], from magnetic bearings in the nonlinear regime [69] and from friction joints [53].

Conclusions
This work presents a method to obtain fast solutions for high-dimensional rotor-foundation systems subjected to nonlinear forces. The basis of the method consist in, firstly, projecting the system (reduced by the CMS) in the linear eigenspace, and then selecting a master mode to enslave all remaining modes. The master mode should be the one with the smallest absolute value of the real part of the eigenvalue (the slowest mode). The approximate manifolds give the relation between the master and the slave modes, and they are obtained from the solution of nonlinear PDEs. With these relations at hand, the equations of motion for the master modes can be solved, and one can obtain the global responses by solving a single pair of equations.
The method was studied in two different systems: a simple and a complex rotor-foundation system. In both cases, the nonlinearity considered comes from fluid-film forces of the bearings. The AIMM was then compared with the responses obtained with the direct numerical integration of the equations. The results show a great capability of reducing the numerical cost and still retaining good accuracy. Therefore, the AIMM is established as a reliable method to perform nonlinear dynamic analysis in rotor-foundation system.