Shooting/continuation based bifurcation analysis of large order nonlinear rotordynamic systems

This study introduces an improved numerical algorithm that is capable of analyzing nonlinear vibrations and bifurcations of general, finite, large-order rotordynamic systems supported on nonlinear bearings. An industrial rotor generally consists of several sections and stages, but numerical shooting/continuation method has been applied to a simple Jeffcott type rotor instead of complex models due to the computational burden of the numerical procedure; it becomes significant when the rotor combined with nonlinear finite bearing models. Here, some mathematical/computational techniques such as a deflation algorithm and the parallel computing are suggested for acceleration along with the conventional treatment of model reduction scheme. An eight-stage compressor rotor supported by two identical five-pad tilting pad journal bearings (TPJB) is selected as a mechanical model to test the numerical incorporation of the algorithms. The rotor beam is modelled with 35 nodes, 140 DOF based on Euler beam theory, and the fluid reaction forces from the two TPJB are calculated using simplex, triangular type finite meshes on the pads. In the numerical procedure, the shooting/continuation combined with the acceleration schemes identifies the solution curves of periodic responses and determines their stability. The orbital motions of coexistent responses are obtained from the solution manifolds.


Introduction
Most of the rotating machineries in industry are not simple as a "Jeffcott" rotor but are consisted of multiple shafts and multiple discs; therefore, the structural vibrations and gyroscopic effects should be taken into consideration for the motions of the machines. The finite element-based rotor beam model with lumped masses and lumped rotational inertias in finite nodes is typical approach to illustrate the general, flexible, multi-mass rotor shaft for numerical analyses. Nonlinear analysis of a large-order system usually demands lots of computation time and resources to incorporate the large number of state elements in both the rotor-beam and the bearing models. Hence, it becomes necessary to combine with a model condensation technique such as real-/complex-eigenvector modal reduction, MATEC Web of Conferences 211, 18003 (2018) https://doi.org/10.1051/matecconf/201821118003 VETOMAC XIV component mode synthesis, and Guyan reduction, which give computational advantages by choosing important modes out of total degree of freedom (DOF) [1]. Component mode synthesis (CMS) method has been selected one very suitable model reduction scheme for nonlinear system analyses. In the procedure, it remains nonlinear associated dynamic components in physical coordinate, i.e. boundary coordinate, and transfer the other components to internal coordinate, then the system DOF is substantially reduced by choosing important modal elements. In a series of research from [2][3][4], CMS had been proven as an accurate dimension reduction method for transient, nonlinear rotordynamics. Meanwhile, Sundararajan and Noah first attempted to combine CMS with shooting-continuation methods to analyze a large order nonlinear rotor-bearing system [5]. In their study, a six-node shaft beam model supported on plain journal bearings was used as a mechanical model; the model has 24 DOF, and the fluid film pressures on the supports are obtained using the short bearing approximation. They obtained the unbalance response amplitude with respect to spin speeds using continuation and determined the local stability based on Floquet theory. Although they showed a possibility to extend it to large order rotor systems, the used model is simple to represent the modern industrial turbomachinery that normally consists of large number of mass/inertia components and complex bearing supports. In particular, if a finite element-based tilting pad journal bearing (TPJB) is employed to the numerical procedure, it can significantly delay the execution speed. In this present study, mathematical and computational acceleration techniques such as a deflation algorithm and the parallel computing strategy are imported to speed up the conventional shooting/continuation-CMS algorithm. The application of two boosting schemes was introduced and verified the efficiency with a finite TPJB model [6]. The objective of this paper is to extend the solubility of the improved numerical algorithms, i.e. shooting/continuation with deflation and the parallel computing, to industrial turbomachinery applications such as a steam turbine compressor.

Theory
This section introduces the numerical methods to obtain single/multiple solutions, bifurcations, and stability of a large order nonlinear rotor-bearing system. This research employed a fixed-interface CMS to condense the large-order systems. Then, dynamics of the reduced system are investigated using a shooting/continuation algorithm, which enhances the execution time through a deflation and the parallel computing strategies.

System condensation
From finite element formulation, a rotor system may have an equation in a form of, Equation (1) can be divided and re-ordered with boundary (nonlinear) coordinate xb, and internal (static) coordinate, xi where Fb is nonlinear force, Fi is force acting at internal coordinate, e.g. unbalance force. Fns,b and Fns,i are non-symmetric element forces from bearings, seals, or other properties at boundary coordinate and internal coordinate, respectively. In CMS, the motion of internal coordinate is assumed as a combination of normal mode and static mode such that: where q is the modal vector. This treatment can transfer internal coordinate to modal coordinate such that: By selecting important modes, which are generally low natural frequency modes under/slightly above the operation speed, total DOF of the original system can be remarkably reduced. Detailed numerical formulations of CMS are well described in [1,5].

Shooting/continuation method
A shooting method identifies a periodic solution of a two-boundary value problem. If a system has τ-period external excitation that is normally corresponded to the rotor spin period, a solution typically exhibits to a rational multiple of the period. Such a nonautonomous case, the shooting searches the solution through the following iterative step such that: where J is Jacobian matrix and I is identity matrix. For a large order system, an initial state vector, x0, for shooting consists of physical and modal coordinates. In physical coordinate, an initial condition (IC) for position and velocity of the journal center at nonlinear bearing locations and IC for bearing itself, for instance, pad angles, is randomly produced within a predefined limit range. On the other hand, a guess for modal coordinate state vector, is generated in descending order. Using a group of initial guesses, the shooting with CMS begins the solution search procedure. If a single or multiple solution is identified, then each solution can stretch out a solution branch with regard to a control parameter, η, using the arc-length continuation method. The numerical procedure can be expressed as follows: where g is a constraint function, and s is the arc-length of the curve. Numerical details of shooting and continuation are well defined in [7].

Parallel computing and deflation
Though a system has the drastically reduced DOF, a shooting/continuation still has an issue to handle the strong local nonlinear force elements such as bearings, seals, squeeze film dampers. As modern turbomachinery tends to increase the speeds and payloads, the bearing support also has to extend performance, which generally requires functional, geometrical complexity such as a tilting pad journal bearing. In such a case, the use of infinitely short/long bearing approximation that was usually applied to the classical shooting/continuation has a limit to accurately predict the fluid reaction force. Alternatively, a finite element-based bearing model seems a proper substitute, it, however, demands lots of computational resource and time to the iterative search procedure in the numerical methods. A combination of acceleration techniques is introduced as a treatment to this issue as shown in Fig. 1. The improved solution procedure utilizes the parallel computing and a deflation algorithm, so Jacobian matrix is provided in parallelized calculations and updated by a newly found solution. The efficiency of the acceleration techniques is quantitatively examined in [6].

Eight-stage compressor supported on TPJBs
The rotor geometry for the present study follows the eight-stage centrifugal compressor model as shown in Fig. 2 (a), which was introduced in [8]. The compressor rotor is modelled as a 35-station finite element beam; the mechanical cross-sectional parameters are specified in [8]. In this study, the added properties at nodes (i.e. weight, transverse moment of inertia and polar moment of inertia) have three times increased values than the original data to lend emphasis on nonlinear behaviors; according to [7], a heavily loaded rotor has higher possibility of various nonlinear responses such as sub-/super-synchronous, quasiperiodic, chaotic vibrations. The compressor rotor is supported on two identical five-pad TPJBs as shown in Fig. 2 (b) and the TPJB specification is shown in Table 1. In general, a rotor hardly maintains perfectly balanced condition, it is assumed that an external excitation induced by imbalance on the fourth compressor disc (i.e. eimb=1.5Cb at 16th node) is applied to the system.

Bifurcation analysis
The bifurcation diagrams are obtained using the introduced shooting/continuation with the acceleration schemes of the parallel computing, deflation, and CMS. The control parameter of continuation is the rotor spin speed (rad/sec), and it proceeds for each harmonic solution identified by the shooting method; here, three coexistent periodic responses, e.g. 1/2τ, 1τ, 2τ, are obtained as shown in Fig. 3. Floquet multipliers of each solution are utilized to determine local stability and bifurcation events. To illustrate the solution branch, the maximum and minimum values of the non-dimensional vertical displacements of rotor at the bearing position, y/Cb, are plotted for each periodic solution as in Fig. 4. In low spin speeds, the journal maintains a stable 1× synchronous response and it locates near the bearing clearance area due to the static load of the compressor. After undergoing a resonance near 550 rad/s, a critical bifurcation occurs at 1060 rad/s, which causes the 1× synchronous response loses its stability and appearance of an unstable 1/2×subsynchronous response; since all the periodic solutions are unstable, it can be interpreted that an aperiodic or a quasi-periodic response appears in steady motion and acts as an attractor. Then, the two periodic solutions converge at 1550 rad/s, and it turns to a stable 1× synchronous response. After all, they separate again at 1650 rad/s with losing the stability. Meanwhile, the unstable 2× super-synchronous response appears all along the operation condition; though it is unstable over the control parameter ranges, this result implies that a rotor system with TPJB supports can shortly exhibit a super-synchronous response by any condition of disturbance during operations.

Conclusion
To examine nonlinear characteristics of an industrial turbomachinery such as a multi-stage compressor rotor supported on multi-pad tilting pad journal bearings, shooting/continuation method is combined with the acceleration schemes of component mode synthesis, the parallel computing, and the deflation algorithm. Though the system is consisted of 150 DOF finite model with the strong local nonlinear elements, the introduced numerical method reveals its accuracy and solubility by providing multiple solution states and bifurcation scenarios.