Assessment of reduced models for the detection of modal interaction through rotor stator contacts

Interactions through direct contact between blade-tips and outer casings in modern turbomachines require complex formulations and subsequent expensive computational efforts when the classical finite element method is considered. The construction of reduced-order models through component mode synthesis techniques usually improves the computational efficiency and may be used for fast parameter studies yielding a better knowledge of the phenomena of interest. In this highly nonlinear framework, the present study is dedicated to the investigation of the capabilities of fixed- and free-interface reduction strategies to handle accurately such problems through a realistic 2D model and complements former results involving a direct modal projection with respective strong kinematic restrictions. The equations of motion are solved using an explicit time integration scheme together with the Lagrange multiplier method where friction is accounted for. The presented work discusses the notions of both displacement and motion convergences and the possibility to conduct fast parameter studies with the use of relevant reduction bases. It also shows that kinematic restrictions artificially enhance the detection of modal interactions.


Abstract
Interactions through direct contact between blade-tips and outer casings in modern turbomachines require complex formulations and subsequent expensive computational efforts when the classical finite element method is considered. The construction of reduced-order models through component mode synthesis techniques usually improves the computational efficiency and may be used for fast parameter studies yielding a better knowledge of the phenomena of interest.
In this highly nonlinear framework, the present study is dedicated to the investigation of the capabilities of fixed-and free-interface reduction strategies to handle accurately such problems through a realistic twodimensional model and complements former results involving a direct modal projection with respective strong kinematic restrictions.
The equations of motion are solved using an explicit time integration scheme together with the Lagrange multiplier method where friction is accounted for. The presented work discusses the notions of both displacement and motion convergences and the possibility to conduct fast parameter studies with the use of relevant reduction bases. It also shows that kinematic restrictions artificially enhance the detection of modal interactions.

Introduction
In modern turbomachines such as aircraft jet engines, improved energy efficiency is achieved by controlling the clearance between the blade-tips and stationary surrounding casings so that aerodynamic leaks are minimized. This strategy inherently leads to more frequent occurrences of direct contacts between the blades and the casings that may originate nonlinear vibrations and subsequent structural damages. Rubbing has been investigated in several related studies [1,2,3], but only a few have focused on the modal interaction mechanism [4,5,6,7] that stands as the phenomenon of interest in the present work. This very specific kind of interaction peculiar to structures featuring cyclic and axi-symmetry can arise under certain conditions: 1. Both structures must vibrate with modes having the same number of nodal diameters (see Fig. 1).
2. Both structures must vibrate at the eigenfrequency of that mode.
3. The backward rotating modes in the bladed disk must travel at the same absolute speed as the forward 1 rotating mode in the casing.
These conditions have been summarized in [8] in a very simple formula where Ω is the angular velocity of the bladed disk and quantities f c , f bd and k respectively stand for the eigenfrequency of the casing, the eigenfrequency of the bladed disk, and the number of diameters of the associated vibration modes. The solution Ω crit of Equation (2) is a critical speed of the system that may lead to an unfavorable configuration due to the simultaneous nonlinear modal resonance of the contacting components. As 3. A locked motion similar to the sustained motion, in which an even number of blades stays permanently in contact with the casing. This type of interaction was detected when the number of blades can be divided by k.
Due to the complexity and size of such formulations discretized through the usual finite element approach, computational times may be prohibitive. Accordingly, it is proposed in the present study to extend the work introduced in [5] by building reduced-order models through two different component mode synthesis techniques and conduct a comparison analysis concerning the detection of specific interacting motions. Among the variety of available component mode synthesis methods [11,9,12,13], only those allowing for a direct treatment of the contact constraints in the reduced-order space should be considered. In other terms, the methods which define physical displacements as interface degrees-of-freedom (DoFs) are actually eligible for this study. The Craig-Bampton (CB) method [9] and the Craig-Martinez (CM) method [10] are the most popular techniques that meet this criteria. They are here adapted to 2D planar models of the bladed disk and the casing in an explicit time-stepping procedure based on the finite-difference scheme with Lagrange multipliers in order to account for contact constraints [14,15]. In this context, a convergence study, with respect to the number of component modes within the reduction basis, introduces the notions of displacement convergence and motion convergence.
In the first section of this paper the modeling of the bladed disk and the casing is described. Then, the contact algorithm and the computation of reduced order model with the different reduction techniques are detailed. Finally, in the last section, results obtained for each type of reduced order model are exposed: modal convergence of the reduced order models computed with component mode synthesis is assessed, the detection of interaction motion is presented and convergences in terms of displacement and motion are extensively revealed for each component mode synthesis method.

Modeling
The presented model involves a two-dimensional representation of the structures of interest very similar to the ones introduced in [5]. It can be regarded as realistic since both normal contact and friction forces treatment between the casing and the tips of the curved blades are accounted for. Moreover, as a first approach the shaft supporting the bladed disk is assumed to be perfectly rigid. Accordingly, the axis of rotation is fixed. This assumption is commonly accpeted in the linear kinematic framework (small displacements and small strains).

Bladed disk
As depicted in Fig. 2, a two-dimensional bladed disk comprising N = 22 blades, typical of modern fans in civil aircraft turbomachines, is considered. Blade geometry and respective displacement field are discretized with the usual Euler-Bernoulli straight beam finite element which has no shear deformation. The global curvature of the blade is achieved through an angle a i between adjacent finite elements. The inter blade phase angle is β = 2π/N .
In this context, any out-of-plane motion of the structure is neglected. It is assumed that cyclic symmetry of the bladed disk is a key-feature to detect modal interaction and that in-plane motions are sufficient to initiate this phenomenon. The structural connection between adjacent blades, namely the disk part, is implemented with curved beam finite elements. The bladed disk is clamped on its rotating shaft. Modal damping is introduced in a general fashion. The total number of DoFs for the bladed disk is 748.
A few free vibration modes of the bladed disk are described in Fig. 3 which explicitly depicts the notion of nodal diameter denoted by k in the present study and specific to structures featuring axi-and cyclicsymmetry [16].

Casing
The casing is modeled as a ring and is discretized using 40 two-noded curved beam finite elements, as pictured in Fig. 4. A few free-vibration modes are pictured in Fig. 5 on which the nodal diameters are also represented. A polar coordinate system with unit vectors e rc and e θc is assigned to the casing for the displacement field and geometry definitions. The initial location of node i is given by the point (R c ,θ i c ), where R c stands for the radius of the casing and l c the length of the element. The finite element description involves four DoFs per node u c , u c,s , v c and v c,s . This choice is motivated by [17] which showed that such formulation is locking free 2 . Table 1 summarizes the geometrical and mechanical properties of the model. After discretization of the weak formulation of the equations of motion, the general problem to be solved is stated as follows for the bladed disk where the superscript m stands for master surface and for the casing where the superscript s stands for slave surface. In each equation, M, K and D respectively stand for the mass matrix stiffness and damping matrices, and F c , F ext and u for the contact forces, any external load and the solution vector. The discretized contact forces F c arise from the contact conditions, referred to as the Kuhn-Tucker optimality conditions, in the following form where Γ m c stands for the master surface, meaning the bladed disk. Conditions (5) specify that the gap g(x) separating the two contacting components cannot be negative along the common interface Γ m c . This gap function may be defined as The Kuhn-Tucker theorem states that there exists a normal pressure t N , assumed positive and acting on the common contact interface in order to enforce the non-penetration condition. Simultaneously, the friction evolution law has to be accounted for. Its simplest formulation lies in the Coulomb model, that is for which µ is the coefficient of friction, v T , the tangential slip rate and t T , the tangential contact force. At convergence, coupled conditions (5) and (7) give rise to the construction of F c . In this work, the assumption of high rotational velocities forbids sticking situations between the blade-tips and surrounding casing and greatly simplifies the final solution strategy [5].
This problem is solved using explicit central differences scheme together with an adapted version of the forward increment Lagrangian method [15,14,5]. In this section, u refers indistinctly to u m or u s respectively mentioned in Eq. (3) and Eq. (4).
The solution method is divided into three steps: 1. Prediction of the displacements u n+1 of the ongoing time-step n by ignoring the contact reactions. These predicted displacements u n+1,p are expressed as 2. Determination of the gap function vector g n+1,p between these structures following Eq. (6). A search algorithm identifies all contactor nodes that have penetrated the target domain. Satisfying the impenetrability conditions implies that the final gap functions g n+1 (linearized when necessary) must be positive or vanish, meaning where the corrected displacements u n+1,c are being calculated. C N is the contact constraint matrix in the normal direction. Eq. (9) is rewritten in an equivalent form where g − n+1,p = Proj R Sg − g n+1,p and Sg = size(g n+1,p ). In other words, only the negative terms of g n+1,p , meaning that a penetration has been detected, are kept in g − n+1,p in order to calculate the corresponding corrected displacements u n+1,c .  Correction of the displacements through the calculation of normal contact and friction forces (due to high relative velocities between the casing and the blade-tips, it is assumed that only sliding occurs). It yields the addition of the unknown Lagrange multipliers λ (or identically t N when using this method) in the governing equations such as The new matrix C NT contains the normal and the pure sliding friction constraints. Ultimately, Eqs. (10) and (11) can be recast in such a form that the Lagrange multipliers are solution of Explicit time integration schemes usually require very small time steps specially for stiff problems such as unilateral contact. In addition, the detection of contact for large finite element formulations dramatically slows down the whole process. This legitimates the use of condensation approaches such as modal projections [5] or more elaborated component mode synthesis techniques.

Reduction techniques
In most industrial applications, finite element models involve a high number of DoFs leading to cumbersome computation times. Numerical techniques such as component mode synthesis procedures allow for a dramatic size condensation of the problem to be solved. The latter are usually coupled with substructuring approaches [11] but only the reduction aspect is considered in the present study.
Reduction from the original full size space u of size n to the reduced order space general coordinates q of size m, with m ≪ n can be achieved by selection of any suitable rectangular transformation matrix of size n × m such as: resulting in a much smaller number of generalized displacements. The content of matrix Φ is prone to discussion here regarding the different criteria of interest: (1) the type of interaction detected, (2) computational times (3) treatment of contact constraints directly in the reduced-order space. Equation (3) thus becomesMq withÂ = Φ T AΦ and A ≡ M, D or K.

Modal projection
By assuming that restricting the kinematics of the blades to a limited motion will help in detecting a specific type of interaction [5], displacement vector u is first developed over two k-nodal diameter modes allowing only those travelling waves having k nodal diameters to propagate in the structure. Accordingly, the coordinate transformation matrix reduces to: where φ k1 and φ k2 are two k-nodal diameter free vibration modes. Both the casing and the bladed disk reduced models become 2 × 2 systems of differential equations. Eventhough the very small size of the models implies a significant system condensation, contact forces can not be computed directly in the reduced-order space and require a backward and forward mapping with the finite element space at each time step. By strongly limiting the motion of the two interacting components, it is expected that the detection of an interaction involving k-nodal diameter modeshapes will be possible. Nevertheless, the treatment of the contact forces is not straightforward. The main advantage of the above mentioned component mode synthesis methods lies in their formulation that preserves a few physical DoFs in the reduced order model. The kinematics restrictions offered by the modal projection are lost but the contact treatment is possible directly in the reduced space.

Component mode synthesis
These physical DoFs are chosen so that they match all the DoFs on which contact forces may be applied over time. By doing so the contact algorithm is significantly improved since, at each time step, both the detection of penetrations and the computation of the induced correction on the predicted displacements can be done directly in the reduced space.
Craig-Bampton (CB) and Craig-Martinez (CM) methods are applied considering each structure as independent substructure connected to each other through a highly nonlinear contact interface.
CB and CM methods first require the distinction of the DoFs of the structure within two groups: the internal DoFs and the boundary DoFs. In general, the definition of these groups is closely related to the loadings applied on the structures. As mentioned above, the boundary is chosen so that it contains any DoF that might be supporting contact force during the simulation, i.e. the two DoFs u and v at the tip of each blade for the bladed disk and the u and v DoFs for each node of the casing. As a result, the boundary of the bladed disk contains 44 DoFs. No component mode synthesis method is applied on the casing since all the u and v DoFs of the structure would be kept in the boundary leading to a reduced order model very close from the finite element model in terms of number of DoFs. One may observe that rotational DoF θ is not kept within the boundary. Indeed, this rotational DoF is not required for managing contact. Moreover, component modes associated with θ may be mathematically redundant with component modes associated with u and v and lead to numerical errors or ill-conditioned reduced matrices.
A reorganization of the general problem (3) to be solved is necessary by separating the DoFs in two groups: the internal DoFs (q i ) and the external DoFs (or boundary: q b ). This yields By definition of the boundary DoFs F c,i = 0.
Craig-Bampton method By definition, the transformation matrix Φ of the CB method stores the following modes [9]: 1. Ψ, a truncated set of eigenmodes of the structure where the boundary DoFs are fixed (see Fig.6(a)), commonly named component modes. The number of component modes is controlled via η; 2. Φ L , a full set of static deformations where, for each boundary DoF, a unitary displacement is applied while the other ones are fixed (see Fig.6(b)). Those are usually called constrained modes.
These notations yield: By construction and due to reorganization (16), Eq. (17) can be expanded as follows: Craig-Martinez method This technique generalizes the usual modal projection defined in subsection 3.2.1 by giving access to some physical displacements directly in the reduced space through an approximated participation of the removed high frequency content modes. The theoretical description given in the paper focuses on the impact of the pseudo-static approximation on the expression of the internal DoFs but the Craig-Martinez change of variable may be obtained faster by using explicitly the generalized force vector. Displacements vector u is first projected onto free mode coordinates of the structure by separating low and high frequency content contributions: where Ψ 1 represents the φ first free vibration modes kept in the reduction basis and Ψ 2 , high frequency modes that will be replaced by a static approximation. Let denote where Ω 1 and Ω 2 are the diagonal matrices containing the square of the angular eigenfrequencies associated with modes Ψ 1 and Ψ 2 such as Ψ T KΨ = Ω which leads to Reading Eq. (21), it yields The pseudo-static approximation neglects the dynamical contributionq 2 of the high frequency modes Ψ 2 , in other words Equation (23) is multiplied by Ψ 2 Ω −1 2 and yields which, combined to Equation (22) leads to Finally, by plugging Eq. (25) in Eq. (19), we find or, in a contracted form: where R is called the residual flexibility matrix.
The key feature of CM method is to insert the boundary DoFs in the unknowns of the reduced system by modifying the transformation of Eq. (27). This is achieved by expressing the external forces F as a function of u b of Eq. (16) and q 1 of Eq. (19). Consequently, Eq. (27) becomes The notations is simplified here, R b and R i are blocks of the square matrix R. However, since F only has coordinates on the boundary DoFs (the boundary DoFs are chosen to be the only one being loaded) the equation can be simplified and from now on: R b = R b1 , R i = R i1 and F = F 1 . Equation (28) becomes then The second block of Eq. (30) is This finally leads to the change of variable The component modes of the CM method are the free vibration modes. Some of them are depicted in Fig. 7.
free vibration modes

Results
The bladed disk is rotating and its angular velocity remains constant during each simulation. The initial blade tip/casing distances is 1 mm for all blades. Between time t = 0 s and t = 2 · 10 −4 s a short external k-nodal diameter impulse is applied on the casing in order to initiate contact. Simulations are conducted for k = 2 and k = 3 pictured in Fig. 8 and increasing angular velocities. Friction coefficient is set to µ = 0.2.

Modal projection
Note: in the following, most of the results associated with the detection of interaction motions are pictured with blade-tip/casing distances such as in Fig. 9(a). Each blue plot stands for the distance between the casing and one of the 22 blades of the bladed disk. When parameter k = 2 only 11 plot are visible due to the symmetry of the contact simulation. Contact between the casing and the blade occurs when the blade-tip/casing distance is equal to zero. Consequently, each interaction motion may be easily observed from these figures. The black dashed line stands for the initial clearance between the bladed disk and the casing. All figures are plotted using the same scale in order to ease comparison. Displacements and time have been normalized.
Even though modal projection leads to two-DoFs systems of equations, simulations remain costly because displacements must be projected back to the physical space at each time step to be correctly corrected when penetrations are detected. On the other hand, it offers the opportunity to detect the expected motion for further investigations.
In agreement with previous results [5], three kinds of motions are observed, namely damped ( Fig. 9(a)), sustained ( Fig. 9(c)) and divergent motions. Among the different sustained motions observed, a specific kind named "locked motion" was detected: a few blades come in permanent contact with the casing, as illustrated in Fig. 9(b), and push a forward travelling wave. The latter is connected to the divisibility of the number of blades N by the number of diameters k of the load exerted on the casing [5]. It means that a 3-diameter configuration might lead to locked sustained motions with a geometry of 30 blades while it could only lead to unlocked sustained motions with 29 blades. As a direct consequence, the contact locations will remain at the same blade tips for all times.

Modal convergence
Before carrying out time simulations based on the reduced order model computed with each component mode synthesis method, their convergence toward the finite element model must be checked. In order to do so the focus is made on the eigenfrequencies of the reduced order model and their evolution while the parameters of reduction of each component mode synthesis technique η and φ are increased. The evolution of the error between some eigenfrequencies of the finite element model and the eigenfrequencies of each reduced order model in function of η for the CB method and φ for the CM method are respectively depicted in Fig. 10(a) and 10(b). It is noticeable that the eigenfrequencies of two reduced order model converge toward the ones of the finite

Detection of interaction motions
CB and CM methods are computationally very efficient since contact forces are computed directly with the reduced models. Anyway, the kinematic restrictions on the displacements of the bladed disk are dropped off leading to difficulties in the detection of k-nodal diameter coincidence.
The results presented in this section are limited to the CB method since, for the sake of brevity, the similar results obtained with the CM method are not detailed. The computed reduced order model of the bladed disk implies η = 88 constraint modes. The modal projection of the casing over its two first k-nodal diameter free vibration modes is considered with k = 2 and k = 3. Consequently, the model of the casing is the same as in section 4.1, it is obtained by modal projection and contains two DoFs associated with each free vibration mode used. Similarly to the previous results, three kinds of motion are detected: damped motions such as the one pictured in Fig. 11(a), sustained motions as displayed in Fig. 11(c) and locked motions depicted in Fig. 11(b). As expected, the frequency content of the vibrations is much richer in the current configuration since the reduced space is spanned by more basis vectors. Consequently, during a locked motion, the noncontacting blades feature independent vibratory patterns in Fig. 11(b) in opposition to a constrained modal projection. During a sustained motion, the amplitudes of vibration are larger as shown in Fig. 11(c).
The previous observations may also be analyzed as follows: a vibration of the casing on a k-nodal shape may be a necessary condition for a k-nodal travelling wave coincidence with the facing bladed-disk to occur. One should notice that results presented in this section are dependant on the reduction parameter η as opposed to the results associated with modal projection given in the previous section. Consequently, it is of great interest to assess the convergence of the given results with respect to a variation of parameter η on which focuses the next section.

Convergence in displacement
The primary challenge in constructing accurate reduced order model is the selection of the basis vectors. Capturing the dynamics of a system is usually achieved through a truncated set of specific shapes. Unfortunately, the latter may be inadequate for nonlinear problems since uniform convergence of coordinate transformation (13) for m → n is proved for linear analyses but does not hold for non-linear analyses especially for stiff systems such as the one investigated here. It is therefore mandatory to perform a convergence analysis in an empiric fashion by comparing the results to the direct finite element formulation with an increasing number (η or φ) of basis vectors.
Even though convergence in time and space are connected, we assume here a constant time step h = 2.5·10 −7 s for all simulations. Convergence of reduced order model computed with CB and CM techniques is studied for each of the three interaction motions detailed above with k = 2 and k = 3 (for sustained motions). Reduction parameters η and φ are chosen in such a way that they remain multiples of 22. This condition ensures the consistency of the enrichment of the reduction basis with the number of harmonics of the system 3 .  Convergence is fast in the case of a damped motion as pictured in Fig. 12(a) but one may clearly observe the great sensitivity of the results to the reduction basis in the case of locked and sustained motions in Figs. 12(b) and 12(c). Indeed, the quality of the approximation of blade tip displacements directly depends on the selected reduction basis vectors. As a consequence, the dimension of the reduced-order model will have an influence on the blades locking on the casing and will also impact the vibration levels the other blades. That is the reason why different vibration levels can be observed in Figs. 12(b) and 12(c). However, results obtained with a large reduction basis (η = 220) are comparable with the finite element results.  Results of the convergence of the reduced order model are shown for damped, locked and sustained motions respectively in Fig. 13(a), 13(b) and 13(c). These results are similar to the ones obtained with the CB method. As previously observed, a large number of modes is required to ensure convergence in terms of blade tip displacements in the cases of locked and sustained motions. This is particularly true for the locked motion for which displacement convergence seems slower than with the CB model.

Craig-Martinez method
Results obtained for the two component mode synthesis methods show a slow convergence of the solution with respect to the finite element reference solution with an increasing number of basis vectors.
For the CM method, it is noticeable that for a small φ (φ = 0 and φ = 44) the results are completely different from the reference solution. In comparison, for similar η, the CB method leads to better results.
A lower threshold for accurate results is η = φ = 88 for damped motions, in the case of sustained or locked motions, even η = φ = 220 do not lead to satisfying results. Comparatively to the size of the finite element model, the condensation ratio for η = φ = 88 -about 18 % -required for displacement convergence may appear too large. It is then worthy to note that this ratio actually reflects that only four families of 22 original bladed disk vibration modes are required for displacement convergence. The ratio between the number of modes of each family of original bladed disk vibration modes in 2D is about 3% 4 while the same ratio for typical 3D industrial models approximately is 0.005%. As a consequence, the condensation ratio for 3D industrial models is expected to be considerably smaller and the combination of component mode synthesis methods with the contact algorithm shall be more efficient. However, it should be underlined that 3D models do not providedue to inconceivable computation times -access to a reference solution as mentioned in [18].
The conclusion of this convergence study is the extreme sensitivity of the results toward the choice of the reduction basis. Since a typical displacement convergence requires very high reduction parameters, it may be of great interest to assess the component mode synthesis methods used with another criterion such as the motion convergence introduced in the next section.

Convergence in motion
As underlined in the previous section, a convergence in terms of displacement is not easily reached due to the inherent difficulties brought by the contact formulation. The fact that the structural displacement field is not exactly retrieved because of the reduction process, due to a phase shift of the blades in contact with the casing, does not mean that a more global convergence with respect to the type of motion detected cannot be obtained. This more global convergence is defined as motion convergence and is now inspected. No finite element solution is calculated for reference because of its prohibitive computation times. Iterations along Ω and µ imply about 700 simulations in order to draw a motion map such as the one depicted in Fig. 14 which is extremely timeconsuming. As a consequence, only the case k = 2 is investigated in this section and only one type of interaction motion may occur: the locked motion.
All the motion maps depicted in this section are obtained with identical simulation parameters and the critical speed defined in (2) is the same for every map and is represented with a red line. The only variable is the type of reduction technique and its associated reduction parameter. Results are presented for both modal projection and component mode synthesis methods for which the convergence is asymptotically assessed by increasing η and φ, η = φ = [22, 44, 66, 88]. First, the map pictured in Fig. 14 highlights the interaction motions detected with modal projection. Three distinct areas may be seen: (1) a damped motion area for low rotational velocities, (2) a locked motion area and a (3) divergent motion area. The boundary between the interaction motion area and the divergent motion area is clearly dependant on the friction coefficient µ. On the contrary, the limit between damped motion area and interaction motion area is unclear in Fig. 14. Although, it is noticeable that the critical speed appears as a higher threshold of the damped motion area. From η = 22 to η = 88, one may observe the great similarity of the four motion maps. This similarity is highlighted by the superposition of the boundaries between interaction motion and divergent motion areas in Fig. 17(a). For each of these four motion maps, the critical speed defines quite precisely the boundary between damped motion and interaction motion areas. The stability of the motion map while η increases highlights that motion convergence is easily reached with the CB method and contrasts with the sensitivity of the displacement convergence observed in the previous section. Motion maps obtained with the CM component mode synthesis method are comparable to the ones obtained with the CB method: the boundary between divergent motion and interaction motion areas is very stablewhich appears clearly in Fig. 17(b) -and only few differences are observable around the critical speed, meaning the boundary between damped motion and interaction motion areas. Same as for the CB method, these maps show that motion convergence is easily reached with the CM method.
Both motion maps obtained with the CB and the CM methods underline that modal projection tends to artificially favor the detection of modal interaction [19]. On the contrary, the matching of motion maps obtained with two different component mode synthesis methods -implying two different reduction bases and reduced spaces -lend weight to the capacity of reduced order models to be used for detecting modal interaction.

Concluding remarks
A combination of component mode synthesis methods with a contact algorithm based on the Lagrange multiplier technique has been introduced in this paper. The use of three types of reduction methods -modal projection over the two first k-nodal diameter free vibration modes and two component mode synthesis methods namely the Craig-Bampton and the Craig-Martinez methods -emphasized the consistency of interaction detection since results presented in this paper show that two different component mode synthesis methods lead to similar results.
It was also pointed out that displacements obtained for each DoF of the reduced order models are extremely sensitive to the modal reduction basis considered due to the highly non linear contact case. Consequently, the notion of motion convergence was introduced in addition to the well-known notion of asymptotic convergence in terms of displacements.
For 2D planar finite element models, it was shown that only a few modes in the modal reduction basis are necessary to determine accurately the different interaction motions. However, in order to get a good accuracy in terms of displacements, the size of the modal reduction basis has to be increased. Results show that Craig-  Martinez component mode synthesis method gives better results with a smaller reduction basis. Finally, it was highlighted that strong kinematic restrictions induced by modal projection artificially enhances the detection of modal interaction.
The fact that a combination of the contact algorithm with component mode synthesis methods lead to consistent results in terms of modal interaction detection is promising for 3D contact simulations to come.