Abstract
The discretisation of rotordynamic systems usually results in a high number of coordinates, so the computation of the solution of the equations of motion is very time consuming. An efficient semianalytic time-integration method combined with a substructure technique is given, which accounts for nonsymmetric matrices and local nonlinearities. The partitioning of the equation of motion into two substructures is performed. Symmetric and linear background systems are defined for each substructure. The excitation of the substructure comes from the given excitation force, the nonlinear restoring force, the induced force due to the gyroscopic and circulatory effects of the substructure under consideration and the coupling force of the substructures. The high effort for the analysis with complex numbers, which is necessary for nonsymmetric systems, is omitted. The solution is computed by means of an integral formulation. A suitable approximation for the unknown coordinates, which are involved in the coupling forces, has to be introduced and the integration results in Green's functions of the considered substructures. Modal analysis is performed for each linear and symmetric background system of the substructure. Modal reduction can be easily incorporated and the solution is calculated iteratively. The numerical behaviour of the algorithm is discussed and compared to other approximate methods of nonlinear structural dynamics for a benchmark problem and a representative example.
1. Introduction
The presence of skew-symmetric matrices is typical for rotordynamic systems. Hence the computation of rotordynamic systems differs from ordinary structural problems due to gyroscopic and circulatory terms, represented by skew-symmetric matrices in the equations of motion; see, for example, Krämer [1]. Modal analysis of such nonsymmetric systems usually involves very expensive computations with complex variables. Nonlinear restoring forces in rotordynamic systems frequently come from the bearings. There is a need for numerical routines for treating the case of catastrophic loading of such structures in the time domain. Frequently nonlinear rotordynamic systems therefore are used as benchmark problems for the verification of the behaviour of new computational solution strategies. In order to keep the computational effort lower, a substructure method is given in the present contribution. The total set of equations of motion is substructured in a way that the slave degrees of freedom have no corresponding components in the vector of the nonlinear restoring force. These equations represent substructure 1 and describe the motion of the linear rotor in the state space. The solution is computed using the configuration-space modes and modal impulse response functions of the linear symmetric parts of the system matrices. The master degrees of freedom are characterized by the presence of nonlinear components. This method of partitioning of the equations of motion of the total system with respect to the nonlinearities of the system is known from the Finite Element Method; see Bathe and Gracewski [2], where usually the Newmark or related direct numerical integration schemes are applied.
In contrast a state-space formulation in connection with a Duhamel-type time integration is given in this contribution. Modal analysis is performed in the confi-guration space considering linear symmetric operators, in order to derive a simple closed-form representation of the two transition matrices of the substructures in the state space. This avoids the more costly analysis with complex numbers due to nonsymmetric matrices and relates the state-space formulation to modal analysis with real eigenvectors of the configuration space. The time integration for both substructures involves the Duhamel-type convolution integral, where the transition matrices are used as kernels. The nonlinearities, the circulatory terms, and the coupling forces are treated as induced forces of the linear system with symmetric matrices. The solution of the second substructure is formulated in the same way. The derivation of an algorithm with a minimum number of nonlinear equations relates the present paper to the component mode analysis, which has been developed by Tongue and Dowell [3] with respect to the transient behaviour of systems with localized nonlinearities and has been extended by Dowell [4]. The Lagrange multiplier method is used in both papers to incorporate the compatibility conditions of the substructures. Kim and Park [5] applied this component mode analysis to nonlinear transient vehicle dynamics. In contrast to the present method the formulation in Dowell [4] is set up separately for the subsystems. Furthermore, we end with a nonlinear system of equations for the increments of the master degrees of freedom, which is solved with the modified Newton-Raphson-method. The present algorithm is formulated in the state space, which is conventional in control engineering, and is suitable not only for nonlinear but also for gyroscopic and circulatory terms. When inserting the coupling conditions and approximating the induced forces, a minimal set of nonlinear equations is derived. Using a time-stepping procedure, this system of incremental equations can be solved iteratively. In order to improve the efficiency of the algorithm different time steps can be applied for each substructure, which can result in a more efficient algorithm, which will be analysed elsewhere. Computing linear and nonlinear rotordynamic effects by means of modal analysis of the linear symmetric parts only, the present paper extends Holl and Irschik [6] with respect to different approximations of the coupling and induced forces and less restrictions for the system matrices. The approximation of the actual time evolution of the physical coordinates is used in the convolution integrals. An extension has been developed by Holl [7] for the case of mass coupling of the substructures. In the more simple case of nonlinearities depending only on the coordinates of substructure 2, the state-space representation can be reduced to a formulation in the configuration space. Further reductions of the computational effort can be achieved using modal analysis with a reduced base of eigenvectors, which is analysed in Holl [8]. As a numerical benchmark application a rotordynamic system is subject to a defined horizontal acceleration component at the bearings.
2. Substructuring with Respect to Nonlinear Restoring Forces
A damped gyroscopic nonlinear system with circulatory forces is considered in the configuration space in the following global form, see, for example, Meirovitch [9]:
, , and are symmetric matrices. The global stiffness matrix is only positive semidefinite if there are rigid body degree-of-freedoms present and the global mass matrix is positive definite. The damping matrix is considered to be of the Rayleigh (proportional) type. If there is any nonclassical damping present in the mechanical system, then for the sake of simple notations these corresponding terms are included in the skew-symmetric gyroscopic matrix here. In Holl [10] the effect of nonclassical damping of the system is studied without substructuring.
It is assumed that the nonlinearity is restricted to the second substructure with the index 2, so that we get
With the global configuration vector and the external forces this partitioning leads to the following matrix-formulation of the two formally decoupled subsets of (1):
where and are the coupling forces of the substructures. Equation (3) may be explained physically in the sense of the structural deformation method, where the equations of motion of the substructure formally can be set up by treating the other parts of the structure to be rigid and fixed; see, for example, Meirovitch [9], Gasch and Knothe [11], and Leung [12]. As can be seen from (3) any mass coupling between the two substructures under consideration is not considered here; this extension is demonstrated in Holl [7] and will be analysed in more detail with special emphasis to the numerical characteristics of the resulting algorithm elsewhere. Frequently the application of lumped mass matrices is preferred in order to prevent a coupling of the mass matrices.
Equation (3) is the starting point of a substructure technique with special emphasis on the system nonlinearities. For rotordynamic systems, usually the partitioning of the equations is done, so that the rotor forms the linear substructure 1, while the nonlinear effects are gathered in substructure 2. As can be seen from (3), in the present contribution the restriction that and are located in two different substructures and do not produce any coupling of the two substructures is dropped, which was used in Holl and Irschik [6]. The coupling forces in (3) become
where the state-space notation
is used and the coupling matrices are
The above mentioned conditions and of the dynamic system are considered. The case of an inertial coupling in (3) can be incorporated using integration by parts; see Holl [7], where different approximations with respect to the integral kernel are reported elsewhere.
3. Formulation for the Substructure with Linear Restoring Forces
The present method takes advantage of the well-known solutions of linear vibration problems with symmetric matrices and proportional damping. A reformulation of (3) for substructure 1 leads to
where the system is excited by the external excitation, the coupling forces, and the motion of the system itself. The mechanical behaviour of substructure 1 is characterized by the symmetric matrices , , and . The increment of the state-vector of (5) at discrete time-intervals is found in the form of the integral equation:
The mass matrix is positive definite and the inverse in (8) needs not to be computed, which is demonstrated in Holl and Irschik [6]. denotes the state vector at the beginning of the time interval under consideration. An efficient computation of the transition matrix was defined in Holl and Irschik [6] by = . This result is based on the consideration of an initial value problem = as the computation of the series expansion of is uncomfortable. In order to get the modal transition matrix, the symmetric eigenvalue problem corresponding to the left hand side of (7) has to be solved:
The eigenvectors are real and orthogonal with respect to and and orthonormalized with respect to , see, for example, Meirovitch [9]: and with . is the identity matrix. The modal expansion of the state-vector is defined by
denotes the vector of modal coordinates. The modal transition matrix mentioned above consists of four diagonal submatrices:
where the diagonal elements of these submatrices are
with the modal parameters
The integrals of (8) can be evaluated now if the time-evolution of the forcing terms , , and would be known. Considering a suitable time-step, which is chosen for the computation of the solution of the nonlinear dynamic system in (8), the time-evolution of these unknown portions of the forcing terms, but not the integrands as a whole, has to be approximated within the time interval. Different approximations seem to be suitable up to second order in the local time coordinate within the time step. For the analysis given here the following approximations are used:
Based on these approximations the integral in (8) can be solved analytically. The coupling forces from (4) are inserted before the integration is performed. For the interpolation of the coordinates of both substructures (12) is used. The solution of the convolution integral then results to
The matrices , and are computed for the cited parabolic approximation:
In (14) the approximations are smoothened by solving the convolution integral. Note that only those submatrices of are involved, which correspond to initial impulses. (8) now becomes
where
takes account for the imposed excitation loading. The algebraic manipulations of putting over to the left-hand side of (15) and solving for finally lead to the compact form:
The form of the used dummy matrices depends on the applied approximation of the time evolution of the coordinate, and for the given case of (12) these matrices are computed by
The matrix in (17) results from the gyroscopic and circulatory terms at the right hand side of (7). Usually the time-step is kept constant during the coarse of computation. In this case the dummy matrices occuring in (18) are calculated only once, before the beginning of the time-stepping procedure. In case of a linear substructure with a large number of degrees-of-freedom, the modal expansion of (10) can be performed with a reduced base of eigenvectors. An analysis involving modal synthesis can save considerable calculation time with only a small reduction in computational accuracy, see Holl [13] for the case without substructuring. For the application of the modal synthesis technique to this method further investigations are performed.
4. Formulation for the Substructure with Nonlinear Restoring Forces
The equations of motion of substructure 2 again are represented as a linear system with symmetric matrices which are excited by external, induced, and nonlinear forces:
The increment of the state-vector is found from (17) simply by changing the indices and extension for the nonlinear restoring force:
Finally the increment of the displacement due to the nonlinearity is computed by
where it is assumed that the nonlinear force varies linear within the time step.
Again the dummy matrices are computed from (18), where the indices have to be changed properly. Equation (17) has to be premultiplied by the inverse of . This matrix inversion is necessary in the following, as the increment of the state vector of substructure 1 of (17) is inserted in (20). The dimension of this matrix can be reduced essentially if modal reduction is applied, which again makes the procedure more efficient. Again is computed once before the beginning of the time stepping procedure. After combining the formulations of both substructures the present nonlinear substructure technique results in an incremental algorithm with a minimum number of nonlinear equations. The reformulated (17) in incremental form is inserted into the master (20), which leads to
where the dummy linear system matrix
and the dummy load vector
are known at the beginning of the time step. Additionally the matrix has the dimension and has to be triangularized. Equation (22) is suitable for general nonlinearities of the type and can be solved using any appropriate iteration procedure for a system of nonlinear equations. Here the modified Newton-Raphson method is implemented, see Bathe [14]. The increment of the vector of the nonlinear force in (22) is computed by
where is expressed by (8) and the modified (17) as a linear function of , which is known from (22).
When including an adaptive time stepping procedure, see, for example, Crisfield [15], all the dummy matrices used in the above algorithm have to be computed after a change in the time step, which is similar to the direct numerical integration procedures. It is mentioned that the modal analysis for the substructures has to be computed only once.
5. Stability Analysis of the Procedure
The analysis of the numerical behaviour of this substructure method involves the numerical dissipation, numerical dispersion, and stability, see, for example, Hughes [16] and Holl [13]. In the following the stability characteristics of the method are demonstrated. Equations (17) and (20) can be reformulated with the additional assumption that the external excitation and nonlinear forces are zero:
is the resulting approximate transfer matrix. A comparison with the corresponding exact transfer matrix formulation
shows the effect of the approximation of a part of the integral. The stability analysis for this procedure defines that the absolute value of the eigenvalues , which is also called spectral radius , has to be less than one: . This result also follows from the corresponding Jordan form of the transfer matrix, see Bathe [14], when considering (26) for time steps, see Holl [13]. For the present approximations the formulation of (26) can be derived and the eigenvalues of this approximate transfer matrix are computed. For sake of simplicity an example with two degrees of freedom is considered, which is described by the matrices , and . The approximate transfer matrix was computed and analysed with respect to the spectral radius. Figure 1 shows the results for the absolute values of the eigenvalues for different dimensionless time steps. The exact solution of the eigenvalues is marked by dashed lines. The good correlation between the exact and the approximate results can be seen in the considered range even for very large time steps. The parameter is calculated from .
6. Numerical Example
The presented method is applied to the rotordynamic system sketched in Figure 2. The rotor is made of steel, with Young’s modulus , Poisson’s ratio 0.3, and mass density . The stationary angular speed of the rotor is 1500 rotations per minute. Equation (1) is derived using Finite Elements for the rotordynamic system, where stiffness and gyroscopic matrices are taken for rotating beam elements and rotor elements. As mentioned above a lumped-parameter mass matrix is used, see Leung [12]. Assuming short bearings, the nonlinear bearing behaviour is described according to the formulas of Childs [17], where is taken for the steady state motion. This steady state is disturbed by the sudden occurrence of a horizontal acceleration component at the bearings, which is given in Figure 3.
Results of the numerical computations with the above semianalytic procedure are shown in Figures 4, 5 and 6, where the time-evolution of the motion at the left and right disk of the rotordynamic system and at the right bearing is shown. The actual displacements are shown in dimensionless form. They are related to the static displacement due to the weight of the rotor itself. The motion is drawn in a plane perpendicular to the rotor axis and for a few cycles of motion. A time-step of 0.002 millisecond has been used so that a converged result is guaranteed.
A comparison to the unconditionally stable version of the Newmark method shows that the effort with the present substructure algorithm is about 20% less than that with the Newmark method, when converged results are considered.
7. Conclusion
The application of the substructure technique in the presented procedure results in a semi-analytic method, which allows for efficient computation of nonlinear rotordynamic systems. Only a few restrictions to the parameters of the system have been introduced. The system matrices can be very general, as for example, the nonproportional damping is considered in the generalized gyroscopic matrix. The stability of the resulting method is shown for the case of a benchmark system. Modal reduction can be implemented easily in the resulting formulation as the nonlinearity is treated in a suitable manner.
Nomenclature
| : | Global mass matrix |
| : | Global damping matrix |
| : | Global gyroscopic matrix |
| : | Global stiffness matrix |
| : | Global stiffness and circulatory matrix |
| : | Global configuration vector |
| : | Global velocity vector |
| : | Global acceleration vector |
| : | Nonlinear restoring force |
| : | Externally imposed force |
| : | Coupling forces of substructure 1 |
| : | Coupling matrix of substructure 1 |
| : | Generalized load vector of substructure 1 |
| : | Global state space vector of substr. 1 |
| : | Increment of the state space vector |
| : | Initial value of |
| : | Modal state space vector |
| : | Vector of modal coordinates |
| : | Global transition matrix |
| : | Transition matrix of substructure 1 |
| : | Approximate global transition matrix |
| : | Modal transition matrix of substructure 1 |
| : | Transition matrix for step load |
| : | Transition matrix for ramp load |
| : | Transition matrix for parabolic load |
| : | Dummy matrix |
| : | Dummy matrix |
| : | Dummy matrix |
| : | Dummy matrix |
| : | Dummy system matrix |
| : | Total modal matrix of the substructure 1 |
| : | Modal matrix of the substructure 1 |
| : | Identity matrix |
| : | Diagonal matrix of the step response |
| : | Diagonal matrix of the impulse response |
| : | th modal eigenfrequency of substr. 1 |
| : | th modal damping of substructure 1 |
| : | th damped modal eigenfrequency |
| : | Degrees of freedom of substructure 1 |
| : | th eigenvalue of the substructure 1. |
Acknowledgment
Support of the author by the Linz Center in Mechatronics (LCM) and the Austrian Center of Competence in Mechatronics (ACCM) is greatly acknowledged.