Abstract
As the use of digital twins for the purpose of structural health monitoring increases, so does the demand for high-quality models. In the area of rotating machines, this means that the casing must also be taken into account. These are often geometrically complex assemblies that require the use of finite element models. The models are then reduced in size by established reduction methods like the Craig–Bampton (CB) reduction in order to keep the solution time, e.g., for a forced response analysis as low as possible. In case of nonlinear contact forces, e.g., friction in bolted joints, a secondary reduction step has to be applied. Here, three different interface reduction (IR) methods are investigated and used for a rotor–casing assembly including frictional damping. The necessary basics and terms are presented. The performance of the methods is evaluated for a linear model on substructure and assembly level and for a nonlinear model on the assembly level. The implementation of the interface reduction within the harmonic balance method (HBM) is presented and tested for different unbalance excitation cases.
1 Introduction
In the field of rotordynamics, the structure that supports the bearings and the rotor plays an increasingly important role due to the growing use of model-based structural health monitoring [1]. In those monitoring applications, the sensors that pick up the vibrations of the machines are often placed on the casing, and thus the transfer function from the rotating system does include the casing. Another aspect is the influence of the dynamic interaction between the rotor and the stator. In order to infer the condition of the rotor, the whole system, including the potentially nonlinear characteristics of the casing, has to be considered.
Usually, the geometry of casings is quite complex, and therefore the structural matrices result from a finite element discretization containing a high number of degrees-of-freedom (DoFs). An established method to handle the large dimension of the matrices is provided by the component mode synthesis [2,3]. The structure is divided into several substructures which are then each separately reduced by, e.g., the Craig–Bampton (CB) or the MacNeal/Rubin-Martinez method. After the reduction, the substructures are assembled by the remaining coupling DoFs that were retained in the reduction. The fact that the coupling DoFs are excluded from the reduction and thus stay in the physical domain allows the use of physically meaningful contact force models, but at the same time, this is also a drawback of the two methods. In case of a casing structure, the selection of the substructures is naturally given by the casing components, and thus the coupling DoFs for the component mode synthesis are the ones that lie in the interfaces of the components. These interfaces can be geometrically large and highly discretized from which it follows that the size of the reduced order model is mainly given by the number of coupling DoFs. If the contact force model is now selected to be nonlinear, e.g., to take into account the damping properties of bolted joints [4,5], the resulting nonlinear equations of motion are too large to be efficiently solved.
Interface reduction (IR) techniques are being developed (and have been developed mainly for contact surfaces of turbine blades) to address this problem [6–9]. In Sec. 2.3, a selection of interface reduction methods is presented and in Sec. 3 these methods are applied to a generic casing model in combination with an unbalanced elastic rotor. A comparison is carried out for a linear and a nonlinear rotor–casing assembly to be able to determine whether the methods can be sensibly used for such systems.
2 Theory
In this section, the theory behind the reduction concepts is presented starting with a short overview of the CB reduction. Then, the different interface reduction methods are described as well as a nonlinear contact force model to consider frictional damping in the interface. Finally, the solution framework using the harmonic balance method (HBM) is outlined.
2.1 Equation of Motion.
The analysis of structural vibration problems is usually based on the mathematical description of the structure which is given by a system of second-order differential equations. Since the basis for this study is a rotor–casing assembly, a distinction must first be made between these two subsystems.
Casing.
where the structural matrices and of each casing component are provided by a finite element program. The matrices are arranged as block diagonal matrices, and the vector of DoFs containing the time-dependent displacements is partitioned accordingly. Since the aim is to model the damping properties of the casing by frictional contact forces in the interfaces, the components in Eq. (1) are not yet connected, and no linear damping matrices are considered.
Rotor.
where the first two entries correspond to translational and the last two to rotational motion. These four DoFs per node are sufficient for the investigation of bending vibrations caused by unbalance forces. The orientation of the coordinate system can be seen in Fig. 1 and chosen, so that the x- and y-axes point in the radial directions and the z-axis coincides with the rotational axis around which the rotor spins with the rotational speed . In Eq. (2), is part of the gyroscopic matrix without the rotational speed . For a simple disk with the four DoFs introduced in Eq. (3), the matrix is presented in Fig. 1 on the right side. The entries are the moment of inertia around the z-axis. The gyroscopic matrix of deflection depends on the internal damping of the rotor which is not considered in this study. Nevertheless, the matrix is retained in this section and in Sec. 2.5 in order to keep the documentation complete.
Assembly.
The force vector contains the contact forces for the coupling of the casing components. In general, nonlinear bearing forces could also be included in if necessary. In case of a linear elastic coupling, the forces can be replaced by stiffness matrices similar to the bearing stiffness matrix . The force vector consists of the forces resulting from unbalances. The in Eq. (4) indicates matrices that have been reduced in advance or are based on reduced matrices. The steps for the calculation of these matrices are presented in Secs. 2.2 and 2.3.
2.2 Craig–Bampton Reduction.
From Eq. (1), the structural matrices of the casing components can be directly extracted. In the following, these components are referred to as substructures. To decrease the size of the full system, each substructure is reduced by a reduction method named after Craig and Bampton. It was introduced in the mid to late 1960s [2,13] and has been a basis for further research ever since [14–18].
The main partition includes all DoFs that should remain in the physical domain after the reduction. In this case, these are the interface DoFs between the substructures and the coupling DoFs for the rotor bearings. All others are assigned to the follower partition.
The size of the resulting matrices in Eq. (12) is specified by the number of fixed interface modes and the number of the main DoFs , which remain completely in .
2.3 Interface Reduction.

Approximation of the real deflection of the interface by superposition of presumed deflection shapes
Rigid Interface.
The new vector of DoFs per interface has now six DoFs that describe the motion of the virtual node and thus also the interface.
Polynomial Basis.
Interface Modes.
The procedure from Eqs. (20)–(24) has to be repeated for every interface of the substructure and for all other substructures forming the complete system.
where contains the orthonormalized interface modes of the interface 1.
2.4 Contact Model.
and corrected if necessary. In case a periodic solution is being sought, the iterative process from Eqs. (28) and (29) has to be done for a few cycles; in Ref. [27], at least three are recommended.

Two-dimensional friction element evaluated for a periodic motion in the tangential plane and a constant normal force. The time at which the state of the spring is depicted is denoted by .
2.5 Harmonic Balance Method.
Preparation.
The resulting vector is the upper part of in Eq. (4).
Nonlinear System of Equations.
The subscript indicates that the contact forces are described in the interface reduction space. In Fig. 5, the minimization of the residual is shown. Starting with an initial guess for , the alternating frequency/time domain method [36] is used to evaluate the nonlinear force in the time domain, see Eqs. (28) and (29). Since the contact model requires the relative displacement in the physical domain, the IR basis has to be used to calculate . Before the contact forces can be inserted into the residual, they need to be transformed back into the interface reduction space by left multiplication with . The iteration process is stopped when the residual falls under a certain threshold. To minimize the work for the solver, the analytical Jacobian matrix is provided, and the secant predictor is used for the initial guess of in case two suitable solutions have already been found. Additionally, the residual scaling is applied to avoid numerical issues due to different magnitudes [37], using the guess provided by the secant predictor as factors.
3 Application
The numerical investigation of the above presented methods is based on the casing and rotor model shown in Fig. 6. The casing in Fig. 6 has an upper and a lower component that are connected by two mesh conforming interfaces. They are simplified by removing any bolt holes in order to make the application of the polynomials as reduction basis fair. The top half (substructure 1) is fixed to ground at the green marked area. The bearings are located on the bottom half (substructure 2). Figure 6(b) shows the rotor with one rigid disk and ten elastic shaft elements. This fairly simple rotor model is sufficient to accurately predict the bending modes that are relevant for unbalance responses [38]. The bearing nodes are located appropriately to match the casing. An overview of the number of DoFs of the three objects is given in Table 1. The CB reduction clearly decreases the number of DoFs significantly. However, the amount of space required to save the mass and stiffness matrices does not decrease in the same order of magnitude because the number of nonzero entries in the structural matrices remains about the same. As the number of main DoFs increases, there will come a point at which the memory advantage is canceled out, but even for this case, the application of IR is necessary.
Number of DoFs after CB reduction and storage space required in the workspace
No. of DoFs | Substructure 1 | Substructure 2 | Rotor |
---|---|---|---|
Initial | 114,921 | 131,934 | 44 |
CB | 1560 | 1564 | 44 |
Interface 1 | 765 | 765 | — |
Interface 2 | 765 | 765 | — |
Initial—space | 173.8 MB | 199 MB | 12.3 kB |
CB—space | 37.14 MB | 74.7 MB | 12.3 kB |
No. of DoFs | Substructure 1 | Substructure 2 | Rotor |
---|---|---|---|
Initial | 114,921 | 131,934 | 44 |
CB | 1560 | 1564 | 44 |
Interface 1 | 765 | 765 | — |
Interface 2 | 765 | 765 | — |
Initial—space | 173.8 MB | 199 MB | 12.3 kB |
CB—space | 37.14 MB | 74.7 MB | 12.3 kB |
Before the interface reduction is carried out, the linear dynamics for the full interface with bonded contacts are investigated in Fig. 7. In Fig. 7(a), the Campbell diagram up to is shown. There are four resonance crossings for the first engine order that are called critical speeds. The corresponding mode shapes are shown in Figs. 7(b)–7(d). The first and second are mainly dominated by the first bending mode of the rotor and the third by the second bending mode. The fourth shows an interaction between rotor and casing: the casing moves side-to-side and the rotor follows a bending motion in the opposite direction. Due to the construction of the bearing support structure and the fixed boundary condition, the stiffness of the casing is anisotropic. This means that all indicated critical speeds could possibly lead to resonance during an unbalance excitation.

Dynamic properties of the assembly with full interfaces and bonded contacts: (a) Campbell diagram, (b) critical speed 1 and 2, (c) critical speed 3, and (d) critical speed 4
3.1 Linear Model.
The different interface reduction methods are first compared based on the individual substructures and the linearly coupled assembly.
Substructure.
is used to evaluate the reduction. The results are shown in Fig. 8 for seven different reduction settings. The digits in the legend correspond to the number of generalized DoFs per interface side. For the GSI modes, the common reduction basis introduced in Eq. (25) is split 50:50 between the two substructures. The first thing that stands out is that the GSI modes are able to approximate the natural frequencies the best with just a few modes, and that they also show good convergence. An increase in the number of interface modes leads to a significant improvement. In case a structure has rigid body modes, like it is the case for substructure 2, the reduction basis of GSI—12 does not contain any elastic interface modes of that substructure; hence, the big differences become obvious in the results for both systems. As expected, the rigid interface and the lower order Legendre polynomials show similar results. These interface modes are just too stiff in order to predict accurate natural frequencies of the individual structures. The storage space required for the structural matrices is listed in Table 2. Even for the largest number of generalized interface DoFs, the decrease in required disk space to save the matrices is substantial.
Assembly.
The substructures are now coupled via the reduced interfaces. For bonded contacts, the rows and columns of the structural matrices that correspond to the relative DoFs are deleted to prevent relative displacement. As a basis of comparison, the first ten critical speeds are used. The error is again calculated based on Eq. (39) and shown in Fig. 9. The results are similar in that the GSI reduced assembly shows the best approximation of the critical speeds, but the difference between the methods is not as large as on the substructure level. This can be attributed to the low participation of the casing in some of the mode shapes, cf. Figs. 7(b) and 7(c).
The bonded contacts are now replaced by linear elastic stiffnesses in the normal and the two tangential directions. The tangential stiffness values are then reduced to simulate a softening behavior as it is often observed for friction contacts with increasing vibration amplitudes. The decrease in the tangential stiffness values is intentionally large in order to have a visible impact in the critical speeds. Figure 10 shows that the first four critical speeds are not strongly dependent on the tangential stiffness. The trend for the error is to increase with decreasing tangential stiffness, and this effect is even more pronounced for the GSI modes than for the LPs. For the seventh and eighth critical Speed, the LPs with order 2 are now slightly more accurate than the GSI modes. Nevertheless, both methods are able to reproduce the dynamics of the assembly well as both have a maximum deviation of around 1%. Finally, the linear unbalance response is calculated. The unbalance force acts on the disk node, and the response is measured at the front bearing. For simplicity and without a loss of generality, it is assumed that the damping for the calculation of the unbalance responses only comes from the viscous bearing damping. The amplitudes shown in Fig. 11 are those of the relative displacement between the rotor node and the corresponding casing node, normalized by the unbalance force amplitude . Up to the third resonance, the frequency response functions are almost identical, even the rigid interface method. Only from the fourth resonance peak, where the rotor–casing interaction is more pronounced, the rigid interfaces are too stiff to accurately predict the critical speed.
3.2 Nonlinear Model.
Based on the findings from the linear investigation, only the RI, LP 2, and GSI 20 will be used for the nonlinear unbalance response. As resonance crossings are usually of special interest, the nonlinear calculation will be limited to the second one. In order to keep the solution times reasonable, only the first harmonic for the HBM is used, so in Eq. (33). Another simplification that is made is the assumption of a homogeneous normal force distribution within the two interfaces. Four different excitation levels are applied by increasing the unbalance mass to predict the behavior of the reduced interfaces for the energy-dependent nonlinear effect of frictional damping. In Fig. 12, the errors (cf. Eq. (39)) of the amplitude calculation at the second resonance in the x-direction are shown. As expected, the RI method is not very accurate in predicting the right vibration amplitudes, as they are generally too stiff. The LP 2 shows better results especially for higher vibration amplitudes, which also fit in with observations from Fig. 10. The best results in terms of calculating the correct vibration amplitudes are achieved by the GSI modes. The convergence data associated with the calculations in Fig. 12 are given in Table 3. All reductions methods reduce the DoFs of the residual, the total amount of contact force calls, which is equivalent to the iteration steps the solver needed to find a solution, and the solution time. Since the solutions of the methods are all different, cf. Fig. 12, the data shown in Table 3 cannot really be compared with each other and instead should give simply an impression of the general performance. Surprisingly, the application of the interface reduction methods does not decrease the amount of contact function calls by a lot but due to the easier handling of the matrices involved in the solution process, still a huge time saving can be achieved.

Error of the nonlinear response amplitude around the second resonance calculated at 300 frequency points
Convergence behavior of the nonlinear calculation from Fig. 12
Interface | No. of DoFs residual | Contact force calls , , , | Time (h) |
---|---|---|---|
Full | 6336 | 1059, 1221, 1504, 1594 | 6.31 |
RI | 264 | 729, 779, 858, 970 | 0.61 |
LP 2 | 432 | 804, 912, 1116, 1274 | 1.19 |
GSI 20 | 376 | 809, 937, 1144, 1238 | 1.31 |
Interface | No. of DoFs residual | Contact force calls , , , | Time (h) |
---|---|---|---|
Full | 6336 | 1059, 1221, 1504, 1594 | 6.31 |
RI | 264 | 729, 779, 858, 970 | 0.61 |
LP 2 | 432 | 804, 912, 1116, 1274 | 1.19 |
GSI 20 | 376 | 809, 937, 1144, 1238 | 1.31 |
Finally, the friction distribution within one of the interfaces is reviewed. The colorbar in Fig. 13 shows the sticking percentage per vibration cycle evaluated at the maximum vibration amplitude. “1” means that the contact is fully stuck at all times, and therefore no friction occurs. The full calculation shows that most of the friction happens at both ends of the interface where the bearing attachment points are located. The RI method cannot reproduce this distribution which explains the poorer performance in the comparison. Slightly better is the LP 2 method which shows less, but still some friction in the middle of interface and more on two ends. On the bottom is the friction distribution for the GSI 20 case shown. This reduction basis approximates the full interface the best.
4 Conclusion
Based on a rotor–casing assembly consisting of one rotor and two casing components connected via two interfaces, three different interface reduction methods were investigated in this paper. To this end, the theory of the methods was first presented, and their implementation within the harmonic balance method–alternating frequency/time domain method-framework for the calculation of nonlinear frequency response functions was demonstrated. It could be shown that the concept of interface reduction can be usefully applied at substructure and assembly level as well as in the calculation of frictionally damped unbalance responses.
Data Availability Statement
The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.
Nomenclature
- a =
active
- bear =
bearing
- c =
constraint, cosine
- case =
casing
- cb =
Craig–Bampton
- cont =
contact
- dyn =
dynamic
- f =
follower
- fi =
fixed interface
- g =
global (includes rotor and casing)
- gsi =
Gram–Schmidt interface
- h =
harmonic component
- i =
interface
- ir =
interface reduced
- lp =
Legendre polynomial
- m =
main
- p =
predict
- ri =
rigid interface
- rot =
rotor
- s =
sine
- =
substructure
- stat =
static
- t =
tangential
- ub =
unbalance