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 [69]. 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.

The starting point for the casing, which often consists of several components, e.g., inlet or outlet housing, is the homogeneous undamped equation of motion
(1)

where the structural matrices (s)M and (s)K of each casing component (s) are provided by a finite element program. The matrices are arranged as block diagonal matrices, and the vector of DoFs q 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.

For the modeling of the elastic rotor, an approach based on the Timoschenko beam theory is used [1012]. Here, the stiffness and the inertia properties of the rotor shaft are approximated by piecewise defined beam sections. In addition, rigid disk elements can be used to model the inertia properties of rotating parts that have little to no influence on the bending stiffness of the rotor. The homogeneous differential equation for a rotationally symmetric rotor reads
(2)
Similar to Eq. (1), the matrices Mrot and Krot contain the inertia and stiffness properties of the rotor, and the matrix Drot contains external damping. The beam sections are defined by two nodes, and each node of the rotor has four DoFs, so that for the nth node the vector is
(3)

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), G0,rot 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 Nrot=ΩN0,rot depends on the internal damping of the rotor which is not considered in this study. Nevertheless, the matrix Nrot is retained in this section and in Sec. 2.5 in order to keep the documentation complete.

Fig. 1
Orientation of the rotor illustrated by a simple disk and the related gyroscopic matrix
Fig. 1
Orientation of the rotor illustrated by a simple disk and the related gyroscopic matrix
Close modal

Assembly.

The combination of both separate systems is realized by the bearings. Here, parallel linear springs and viscous dampers are used for that purpose as a simplified numerical bearing model. The entire equation of motion for the assembly is then
(4)

The force vector fcont contains the contact forces for the coupling of the casing components. In general, nonlinear bearing forces could also be included in fcont if necessary. In case of a linear elastic coupling, the forces can be replaced by stiffness matrices similar to the bearing stiffness matrix Kbear. The force vector fub 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 [1418].

The idea of the method is to split (s)q of the substructure (s) into two groups named main (m) and follower (f) DoFs
(5)
and to replace the follower part (s)qf by a static and dynamic component
(6)

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 first part of Eq. (6) results from the second row of Eq. (5) for the static case
(7)
The matrix (s)Ψc contains the so-called constraint modes which ensure the static completeness of the main DoFs. For the calculation of the dynamic part of Eq. (6), first an eigenvalue problem is set up by the follower partition of Eq. (5)
(8)
(9)
The solution of Eq. (8) is only calculated for a predefined frequency range, and the resulting eigenvectors are collected in the modal matrix (s)Φfi. They are called fixed interface modes because the deletion of the main partition for the setup of Eq. (8) is identical to constraining the main DoFs to be equal to zero. Using the truncated modal matrix (s)Φfi, the dynamic part (s)qf,dyn can be substituted
(10)
by the new set of modal DoFs (s)ηfi which has significantly less entries than (s)qf,dyn. Equations (7) and (10) are finally combined to build the reduction basis in matrix form
(11)
which then is used to reduce the substructure by inserting Eq. (11) into Eq. (5) and a left multiplication with (s)TcbT
(12)

The size of the resulting matrices in Eq. (12) is specified by the number of fixed interface modes Nfi and the number of the main DoFs Nm, which remain completely in (s)qcb.

2.3 Interface Reduction.

If NmNfi, the reduction using (s)Tcb from Eq. (12) is not sufficient, and a secondary reduction step has to be carried out, which is called IR. The concept of IR is similar to Eq. (10) where a number of assumed deflection shapes, in this case solutions of an eigenvalue problem of the follower part of the structure, are used to approximate the actual deflection. The same idea can be applied to the main part of the substructure [19,20] as visualized in Fig. 2. The deflections of the main nodes are composed from an interface reduction basis (s)Tir which decreases the number of physical interface DoFs to a smaller set of generalized interface DoFs (s)ηi. The application to the substructure is identical to Eq. (12). The secondary reduction matrix for a substructure having one set of active nodes, i.e., DoFs that should remain in the physical domain, and two interfaces is
(13)
Fig. 2
Approximation of the real deflection of the interface by superposition of presumed deflection shapes
Fig. 2
Approximation of the real deflection of the interface by superposition of presumed deflection shapes
Close modal

In the past, several methods for the calculation of reduction bases (s)Tir have been developed [68,1923], three of which are presented and compared in this study.

Rigid Interface.

The first method investigated is introduced in Ref. [23]. The authors allow the interface to move as a rigid body, which leaves a maximum of six DoFs per interface side. The motivation for the development of this concept is to efficiently connect a flexible and a stiff body. The calculation of the reduction basis requires to define a virtual (v) node which is then used for the coupling, thus having three translational and three rotational DoFs. In this study, the location of this node is always chosen to be at the geometric center of the interface nodes. Using this location, the distances to the nth node are
(14)
The nodal deflections of the rigid interfaces (RI) can then be expressed by the linearized relations
(15)
Casting Eq. (15) in matrix form, the reduction basis is compactly written as
(16)
where the matrices (s)An contain the kinematic constraints for the rotational DoFs
(17)

The new vector of DoFs per interface has now six DoFs (s)qv=[(s)xv,(s)yv,(s)zv,(s)φx,v,(s)φy,v,(s)φz,v]T that describe the motion of the virtual node and thus also the interface.

Polynomial Basis.

Another method is using orthogonal polynomials as basis functions presented in Refs. [22] and [24]. According to the authors, it can be seen as further development of the rigid interface method since the polynomials include the case of rigid body motion. In Ref. [22], it is proposed to use Legendre polynomials (LP) as basis functions because of their natural orthogonality in the range of [1,1]. This prevents redundancy in the reduction basis which would lead to an unnecessary large number of needed functions to achieve a certain, previously specified limit of convergence. Another positive characteristic is the extension of the polynomials to two-dimensional (2D) surfaces by a dyadic product, shown in Fig. 3. The polynomials can therefore be easily applied to 2D surfaces. Figure 3 also clearly shows that the setup of the polynomials is done on a [1,1]×[1,1] grid. In order to apply the method to a real interface, the polynomials have to be mapped to the actual geometry of the interface. This mapping cannot be done to arbitrarily complex interfaces, which marks a limitation of this method. The reduction matrix is then built from the function values at the node locations oriented as column vectors P, exactly like the matrix that is drawn in Fig. 2. The degree of reduction can be varied for the different spatial directions, e.g., depending on the normal or the tangential direction in relation to the interface
(18)
Fig. 3
Construction of 2D Legendre polynomials based on two one-dimensional polynomials
Fig. 3
Construction of 2D Legendre polynomials based on two one-dimensional polynomials
Close modal

Interface Modes.

Unlike the first two methods, which use mathematically motivated basis functions, the third one makes use of the actual structural information. Based on the characteristic constraint modes introduced in Ref. [20], the Gram–Schmidt interface (GSI) modes were developed and utilized in Refs. [7], [25], and [26]. In order to calculate these interface modes, an eigenvalue problem using the main partition of Eq. (12) is set up. To illustrate the procedure, the same composition of the main partition from Eq. (13)
(19)
is used. That means the substructure (s) has two interfaces and one set of active nodes. The first step in the setup of the eigenvalue problem for the interface number 1 is the rearrangement of the groups in Eq. (19)
(20)
in the way that (s)qm,i1 is the first set in the vector. The matrices corresponding to the main part are sorted accordingly, and the solution of the first eigenvalue problem
(21)
provides the interface modes (s)φi,1,j that are collected in a modal matrix (s)Φi1. The content of the modal matrix (s)Φi1
(22)
not only includes information about the interface 1 but also about the other two main groups, which are merged as indicated by the index r. From Eq. (22), it becomes clear that the main groups are coupled. Rearranging Eq. (22) to eliminate (s)ηr from the expression to calculate (s)qm,i1 leads to
(23)
To achieve the DoFs reduction, the modal matrix (s)Φi1 needs to be truncated, which deletes the marked coupling terms in Eq. (23) and simplifies it to
(24)

The procedure from Eqs. (20)(24) has to be repeated for every interface of the substructure (s) and for all other substructures forming the complete system.

If the two substructures (s) and () are to be joined via a reduced interface, the reduction bases have to be compatible. Since they originate from different eigenvalue problems (cf. Eq. (21)), a common basis has to be built
(25)
from both solutions. The last step is to apply the name-giving Gram–Schmidt orthonormalization to Eq. (25) because the mode shapes that form i1Φ are only a part of the full mode shapes of Eq. (21) and thus can be very similar and far from orthogonal which could lead to numerical issues. The reduction of the main DoFs for interface 1 is then
(26)
(27)

where Tgsi,i1 contains the orthonormalized interface modes of the interface 1.

2.4 Contact Model.

For the investigation of the performance of the IR methods in presence of nonlinear contact forces, the node-to-node 3D-Jenkins model with a constant normal force is utilized [2731]. The presumably complex dynamic behavior of the casing makes it impossible to predict preferred directions in the tangential plane of the interface, so that the coupling of tangential DoFs should be taken into account [32]. In Fig. 4, the concept of the element is presented. It consists of the tangential stiffness kt that is attached to a Coulomb-slider. As long as the spring force is smaller than the friction force fr, the slider stays at the same point. If the friction force is overcome, the slider starts moving into the direction of the applied force, as indicated with the long arrow at the normalized time τ=τp. On the right side of Fig. 4, the tangential forces over one period are shown. It is clearly visible from the amplitude of both forces that the maximum force is limited by the friction force fr. The evaluation of the tangential forces must be carried out iteratively in a predictor–corrector fashion, e.g., the change of the force at time ti is assumed to belong
(28)
where Δt is the relative tangential displacement. This linear elastic prediction has to be checked
(29)

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.

Fig. 4
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 τp.
Fig. 4
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 τp.
Close modal

2.5 Harmonic Balance Method.

The calculation of the nonlinear unbalance responses is carried out using the HBM. It is commonly used in the field of blade dynamics when nonlinear forces due to friction models are involved [33,34].

Preparation.

An assembly consisting of two substructures connected by two interfaces is used for explanation of the solution procedure. The reduced DoFs are
(30)
where the generalized coordinates of the interfaces are labeled with main and follower, just to be able to tell them apart. The structural matrices of the substructures are then assembled to get the global vector of DoFs
(31)
It is now beneficial to replace the follower interface DoFs by relative DoFs Δηi=ηi,fηi,m using a matrix containing only identity matrices
(32)

The resulting vector qrel is the upper part of qg in Eq. (4).

Nonlinear System of Equations.

As part of the application of the HBM to Eq. (4), a periodic oscillation of the time-dependent coordinates qg and contact forces fcont are assumed
(33)
(34)
The number of considered frequency components is denoted with H, and Q and F contain the amplitudes of the sine and cosine parts. The excitation force due to unbalance acting on the nth rotor-node is formulated in the same way
(35)
where mu is an unbalance mass that is located on a radius ru at an angle γ. Separating the sine and cosine parts gives the equations of motion in the frequency domain
(36)
The first part of Eq. (36) can be summarized as dynamic stiffness matrix
(37)
that depends on the rotational speed. The unknowns in Eq. (37) are the interdependent Fourier coefficients Qg and Fcont which are iteratively determined using a nonlinear solver. A residual must be set up for this purpose where only the DoFs with a nonlinear force are included. This can be achieved by a dynamic reduction as suggested in Ref. [35] to the relative DoFs Δηi (cf. Eq. (32))
(38)

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 Qi, 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 Δti. 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 TirT. 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 Qi 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.

Fig. 5
Scheme of the alternating frequency/time domain method including interface reduction
Fig. 5
Scheme of the alternating frequency/time domain method including interface reduction
Close modal

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.

Fig. 6
Model of the rotor–casing assembly used for the numerical studies
Fig. 6
Model of the rotor–casing assembly used for the numerical studies
Close modal
Table 1

Number of DoFs after CB reduction and storage space required in the workspace

No. of DoFsSubstructure 1Substructure 2Rotor
Initial114,921131,93444
CB1560156444
Interface 1765765
Interface 2765765
Initial—space173.8 MB199 MB12.3 kB
CB—space37.14 MB74.7 MB12.3 kB
No. of DoFsSubstructure 1Substructure 2Rotor
Initial114,921131,93444
CB1560156444
Interface 1765765
Interface 2765765
Initial—space173.8 MB199 MB12.3 kB
CB—space37.14 MB74.7 MB12.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 15×103min1 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.

Fig. 7
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
Fig. 7
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
Close modal

3.1 Linear Model.

The different interface reduction methods are first compared based on the individual substructures and the linearly coupled assembly.

Substructure.

Starting on the substructure level, the first ten natural frequencies corresponding to elastic mode shapes are used for the comparison. Since the performance of the IR is to be assessed, the reference is always the CB reduced model having the full interface. The error e of the approximation
(39)

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.

Fig. 8
Influence of interface reduction methods on the natural frequencies on substructure level
Fig. 8
Influence of interface reduction methods on the natural frequencies on substructure level
Close modal
Table 2

Storage required to save the reduced matrices

No. of DoFsSubstructure 1 (kB)Substructure 2 (kB)
LP 0—320.2425
RI—627.5633.6
GSI—1245.5652.56
GSI—2076.5685.56
LP 2—27110.24121
No. of DoFsSubstructure 1 (kB)Substructure 2 (kB)
LP 0—320.2425
RI—627.5633.6
GSI—1245.5652.56
GSI—2076.5685.56
LP 2—27110.24121

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 Δηi 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).

Fig. 9
Influence of interface reduction methods on the critical speeds for bonded contacts
Fig. 9
Influence of interface reduction methods on the critical speeds for bonded contacts
Close modal

The bonded contacts are now replaced by linear elastic stiffnesses in the normal and the two tangential directions. The tangential stiffness values kt 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 muruΩ2. 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.

Fig. 10
Influence of interface reduction methods on the critical speeds for linear elastic contacts
Fig. 10
Influence of interface reduction methods on the critical speeds for linear elastic contacts
Close modal
Fig. 11
Influence of interface reduction methods on the unbalance response function
Fig. 11
Influence of interface reduction methods on the unbalance response function
Close modal

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 H=1 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 mu 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.

Fig. 12
Error of the nonlinear response amplitude around the second resonance calculated at 300 frequency points
Fig. 12
Error of the nonlinear response amplitude around the second resonance calculated at 300 frequency points
Close modal
Table 3

Convergence behavior of the nonlinear calculation from Fig. 12 

InterfaceNo. of DoFs residualContact force calls 1mu, 2mu, 4mu, 6muTime (h)
Full63361059, 1221, 1504, 15946.31
RI264729, 779, 858, 9700.61
LP 2432804, 912, 1116, 12741.19
GSI 20376809, 937, 1144, 12381.31
InterfaceNo. of DoFs residualContact force calls 1mu, 2mu, 4mu, 6muTime (h)
Full63361059, 1221, 1504, 15946.31
RI264729, 779, 858, 9700.61
LP 2432804, 912, 1116, 12741.19
GSI 20376809, 937, 1144, 12381.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.

Fig. 13
Stick–slip distribution within the interface
Fig. 13
Stick–slip distribution within the interface
Close modal

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

Latin Letters
D =

damping matrix

G =

gyroscopic matrix

K =

stiffness matrix

M =

mass matrix

N =

gyroscopic matrix of deflection

q =

generalized degrees-of-freedom

T =

transformation or reduction matrix

Greek Symbols
Δ =

relative value

η =

generalized coordinate of a deflection/mode shape

φ =

eigenvector

Φ =

modal matrix

ω =

natural frequency

Ω =

rotational speed

Superscripts and Subscripts
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

(s) =

substructure

stat =

static

t =

tangential

ub =

unbalance

References

1.
Thelen
,
A.
,
Zhang
,
X.
,
Fink
,
O.
,
Lu
,
Y.
,
Ghosh
,
S.
,
Youn
,
B. D.
,
Todd
,
M. D.
,
Mahadevan
,
S.
,
Hu
,
C.
, and
Hu
,
Z.
,
2022
, “
A Comprehensive Review of Digital Twin—Part 1: Modeling and Twinning Enabling Technologies
,”
Struct. Multidiscip. Optim.
,
65
(
12
), Paper No. 354.10.1007/s00158-022-03425-4
2.
Hurty
,
W. C.
,
1965
, “
Dynamic Analysis of Structural Systems Using Component Modes
,”
AIAA J.
,
3
(
4
), pp.
678
685
.10.2514/3.2947
3.
Gruber
,
F. M.
, and
Rixen
,
D. J.
,
2016
, “
Evaluation of Substructure Reduction Techniques With Fixed and Free Interfaces
,”
Strojniški Vestn. - J. Mech. Eng.
,
62
(
7–8
), pp.
452
462
.10.5545/sv-jme.2016.3735
4.
Böswald
,
M.
,
Link
,
M.
,
Meyer
,
S.
, and
Weiland
,
M.
,
2002
, “
Investigations on the Non-Linear Behaviour of a Cylindrical Bolted Casing Joint Using High Level Base Excitation Tests
,”
Proceedings of ISMA
, Leuven, Belgium, Sept. 16–18, pp.
1203
1212
.https://www.researchgate.net/publication/331167738_Investigations_on_the_Non-xcitation_Tests
5.
Böswald
,
M.
,
Link
,
M.
, and
Schedlinski
,
C.
,
2005
, “
Computational Model Updating and Validation of Aero-Engine Finite Element Models Based on Vibration Test Data
,”
Proceedings of the International Forum on Aeroelasticity and Structural Dynamics
(
IFASD
), Munich, Germany, June 28–July 1.https://www.researchgate.net/publication/267836643_COMPUTATIONAL_MODEL_UPDATING_AND_VALIDATION_OF_AEROENGINE_FINITE_ELEMENT_MODELS_BASED_ON_VIBRATION_TEST_DATA
6.
Tran
,
D.-M.
,
2009
, “
Component Mode Synthesis Methods Using Partial Interface Modes: Application to Tuned and Mistuned Structures With Cyclic Symmetry
,”
Comput. Struct.
,
87
(
17–18
), pp.
1141
1153
.10.1016/j.compstruc.2009.04.009
7.
Battiato
,
G.
,
2017
, “
Vibrations Prediction and Measurement of Multi-Stage Bladed Disks With Non Linear Behavior Due to Friction Contacts
,”
Ph.D. thesis
,
Graduate School of Politecnico di Torino
, Turin, Italy.10.6092/POLITO/PORTO/2680969
8.
Krattiger
,
D.
,
Wu
,
L.
,
Zacharczuk
,
M.
,
Buck
,
M.
,
Kuether
,
R. J.
,
Allen
,
M. S.
,
Tiso
,
P.
, and
Brake
,
M. R. W.
,
2019
, “
Interface Reduction for Hurty/Craig-Bampton Substructured Models: Review and Improvements
,”
Mech. Syst. Signal Process.
,
114
, pp.
579
603
.10.1016/j.ymssp.2018.05.031
9.
Förster
,
A.
, and
Panning-von Scheidt
,
L.
,
2022
, “
Interface Reduction in an Equivalent Linearization Algorithm for Nonlinearly Coupled Systems Under Random Excitation
,”
ASME
Paper No. GT2022-78367.10.1115/GT2022-78367
10.
Nelson
,
H. D.
,
1980
, “
A Finite Rotating Shaft Element Using Timoshenko Beam Theory
,”
ASME J. Mech. Des.
,
102
(
4
), pp.
793
803
.10.1115/1.3254824
11.
Hutchinson
,
J. R.
,
2001
, “
Shear Coefficients for Timoshenko Beam Theory
,”
ASME J. Appl. Mech.
,
68
(
1
), pp.
87
92
.10.1115/1.1349417
12.
Friswell
,
M. I.
,
Penny
,
J. E. T.
,
Garvey
,
S. D.
, and
Lees
,
A. W.
,
2009
,
Dynamics of Rotating Machines
, 1st ed.,
Cambridge University Press
, Cambridge, UK.
13.
Craig
,
R. R.
, and
Bampton
,
M. C. C.
,
1968
, “
Coupling of Substructures for Dynamic Analyses
,”
AIAA J.
,
6
(
7
), pp.
1313
1319
.10.2514/3.4741
14.
Kammer
,
D. C.
, and
Baker
,
M.
,
1987
, “
Comparison of the Craig-Bampton and Residual Flexibility Methods of Substructure Representation
,”
J. Aircr.
,
24
(
4
), pp.
262
267
.10.2514/3.45435
15.
Craig
,
R.
, Jr.
,
2000
, “
Coupling of Substructures for Dynamic Analyses—An Overview
,”
AIAA
Paper No. 2000-1573.10.2514/6.2000-1573
16.
Rixen
,
D. J.
,
2011
, “
Interface Reduction in the Dual Craig-Bampton Method Based on Dual Interface Modes
,”
Linking Models and Experiments
, Vol. 2,
Springer
,
New York
, pp.
311
328
.
17.
Kim
,
J.-G.
, and
Lee
,
P.-S.
,
2015
, “
An Enhanced Craig-Bampton Method
,”
Int. J. Numer. Methods Eng.
,
103
(
2
), pp.
79
93
.10.1002/nme.4880
18.
Gruber
,
F. M.
, and
Rixen
,
D.
,
2018
, “
Comparison of Craig-Bampton Approaches for Systems With Arbitrary Viscous Damping in Dynamic Substructuring
,”
Dynamics of Coupled Structures
, Vol.
4
,
Springer International Publishing
, Cham, Switzerland, pp.
35
49
.
19.
Balmes
,
E.
,
1996
,
Use of Generalized Interface Degrees of Freedom in Component Mode Synthesis
,
IMAC
, Dearborn, MI.
20.
Castanier
,
M. P.
,
Tan
,
Y.-C.
, and
Pierre
,
C.
,
2001
, “
Characteristic Constraint Modes for Component Mode Synthesis
,”
AIAA J.
,
39
(
6
), pp.
1182
1187
.10.2514/2.1433
21.
Herrmann
,
J.
,
Maess
,
M.
, and
Gaul
,
L.
,
2010
, “
Substructuring Including Interface Reduction for the Efficient Vibro-Acoustic Simulation of Fluid-Filled Piping Systems
,”
Mech. Syst. Signal Process.
,
24
(
1
), pp.
153
163
.10.1016/j.ymssp.2009.05.003
22.
Holzwarth
,
P.
, and
Eberhard
,
P.
,
2015
, “
Interface Reduction for CMS Methods and Alternative Model Order Reduction
,”
IFAC-PapersOnLine
,
48
(
1
), pp.
254
259
.10.1016/j.ifacol.2015.05.005
23.
Lindberg
,
E.
,
Hörlin
,
N.-E.
, and
Göransson
,
P.
,
2013
, “
Component Mode Synthesis Using Undeformed Interface Coupling Modes to Connect Soft and Stiff Substructures
,”
Shock Vib.
,
20
(
1
), pp.
157
170
.10.1155/2013/262354
24.
Carassale
,
L.
, and
Maurici
,
M.
,
2018
, “
Interface Reduction in Craig–Bampton Component Mode Synthesis by Orthogonal Polynomial Series
,”
ASME J. Eng. Gas Turbines Power
,
140
(
5
), p.
052504
.10.1115/1.4038154
25.
Battiato
,
G.
,
Firrone
,
C. M.
,
Berruti
,
T. M.
, and
Epureanu
,
B. I.
,
2018
, “
Reduction and Coupling of Substructures Via Gram–Schmidt Interface Modes
,”
Comput. Methods Appl. Mech. Eng.
,
336
, pp.
187
212
.10.1016/j.cma.2018.03.001
26.
Battiato
,
G.
, and
Firrone
,
C. M.
,
2019
, “
Reduced Order Modeling of Large Contact Interfaces to Calculate the Non-Linear Response of Frictionally Damped Structures
,”
Procedia Struct. Integr.
,
24
, pp.
837
851
.10.1016/j.prostr.2020.02.074
27.
Sanliturk
,
K. Y.
, and
Ewins
,
D. J.
,
1996
, “
Modelling Two-Dimensional Friction Contact and Its Application Using Harmonic Balance Method
,”
J. Sound Vib.
,
193
(
2
), pp.
511
523
.10.1006/jsvi.1996.0299
28.
Yang
,
B.
, and
Menq
,
C.
,
1998
, “
Characterization of 3D Contact Kinematics and Prediction of Resonant Response of Structures Having 3D Frictional Constraint
,”
J. Sound Vib.
,
217
(
5
), pp.
909
925
.10.1006/jsvi.1998.1802
29.
Genzo
,
A.
,
Sextro
,
W.
, and
Panning-von Scheidt
,
L.
,
2006
, “
Dynamic Behaviour of Elastic Bodies Coupled by Extended Friction Contacts
,”
Proceedings of ISMA 2006: International Conference on Noise and Vibration Engineering
, Vol.
3
, Leuven, Belgium, Sept. 18–20, pp.
1303
1317
.https://www.fis.uni-hannover.de/portal/en/publications/dynamic-behaviour-of-elastic-bodies-coupled-by-extended-friction-contacts(d5a50bec-6686-45b4-9c7c-8c6b8c175e64).html
30.
Sextro
,
W.
,
2007
,
Dynamical Contact Problems With Friction
, 1st ed.,
Springer
,
Berlin, Heidelberg, Germany
.
31.
Bograd
,
S.
,
Reuss
,
P.
,
Schmidt
,
A.
,
Gaul
,
L.
, and
Mayer
,
M.
,
2011
, “
Modeling the Dynamics of Mechanical Joints
,”
Mech. Syst. Signal Process.
,
25
(
8
), pp.
2801
2826
.10.1016/j.ymssp.2011.01.010
32.
Gross
,
J.
,
Armand
,
J.
,
Lacayo
,
R. M.
,
Reuss
,
P.
,
Salles
,
L.
,
Schwingshackl
,
C. W.
,
Brake
,
M. R. W.
, and
Kuether
,
R. J.
,
2016
, “
A Numerical Round Robin for the Prediction of the Dynamics of Jointed Structures
,”
Dynamics of Coupled Structures
, Vol.
4
,
Springer International Publishing
, Cham, Switzerland, pp.
195
211
.
33.
Menq
,
C.-H.
, and
Griffin
,
J. H.
,
1985
, “
A Comparison of Transient and Steady State Finite Element Analyses of the Forced Response of a Frictionally Damped Beam
,”
ASME J. Vib. Acoust.
,
107
(
1
), pp.
19
25
.10.1115/1.3274709
34.
Krack
,
M.
,
Salles
,
L.
, and
Thouverez
,
F.
,
2017
, “
Vibration Prediction of Bladed Disks Coupled by Friction Joints
,”
Arch. Comput. Methods Eng.
,
24
(
3
), pp.
589
636
.10.1007/s11831-016-9183-2
35.
Poudou
,
O.
, and
Pierre
,
C.
,
2003
, “
Hybrid Frequency-Time Domain Methods for the Analysis of Complex Structural Systems With Dry Friction Damping
,”
AIAA
Paper No. 2003-1411.10.2514/6.2003-1411
36.
Cameron
,
T. M.
, and
Griffin
,
J. H.
,
1989
, “
An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems
,”
ASME J. Appl. Mech.
,
56
(
1
), pp.
149
154
.10.1115/1.3176036
37.
Krack
,
M.
, and
Gross
,
J.
,
2019
,
Harmonic Balance for Nonlinear Vibration Problems
, 1st ed.,
Springer International Publishing
, Cham, Switzerland.
38.
Genta
,
G.
,
Silvagni
,
M.
, and
Qingwen
,
C.
,
2013
, “
Dynamic Analysis of Rotors: Comparison Between the Simplified One-Dimensional Results and Those Obtained Through 3-D Modeling
,”
XXI Congresso nazionale AIMETA
, Turin, Italy, Sep. 17–20, pp.
1
10
.https://core.ac.uk/reader/19751430