Modal Derivatives for Efficient Vibration Prediction of Geometrically Nonlinear Structures with Friction Contact

: This paper evaluates the performance of the Rubin reduction methods, enhanced with static modal derivatives, for vibration analysis of geometrically nonlinear structures with friction contact. Static modal derivatives are computed numerically based on Rubin reduction, which includes free interface normal modes and residual flexibility attachment modes, by introducing a finite displacement around these modes. Then, the most relevant static modal derivatives are selected using an improved strategy that incorporates weighting factors derived from both a nonlinear static analysis and a geometrically linear transient run. This enhanced Rubin method is also compared with the previously used enhanced Craig–Bampton method, which is based on fixed normal modes, constraint modes, and their static derivatives. The effectiveness of these methods is demonstrated through vibration analysis of a geometrically nonlinear beam in different contact configurations. Both methods proved their robustness, achieving accurate results with a relatively small number of modes in the reduced space, thus ensuring low online computation times. Furthermore, the analyses show the significant impact of using a geometrically nonlinear model on the accurate prediction of a contact state.


Introduction
The dynamic behavior of slender bodies as aero-engine blades can be characterized simultaneously by large deformation and rubbing phenomena.A typical example is the fan blade of engines with ultra-high by-pass ratios [1].On one side, large deformation activates distributed geometric nonlinearities [2,3], while rubbing phenomena introduce localized nonlinearity to the system [4].These nonlinearities can give rise to complex dynamic behavior, which must be accurately captured during the design process.The Finite Element (FE) method is typically used as an accurate and robust tool for simulating the structural response.However, the nonlinear dynamic analysis of a high-fidelity FE model with a large number of degrees of freedom can be prohibitively expensive.Model Order Reduction (MOR) techniques are employed to reduce this computing cost.In the literature, numerous reduction techniques have been developed to address either geometric or contact nonlinearity separately.Given the distinct nature of these two nonlinearities, the developed reduction techniques are essentially different.This paper focuses on a reduced model capable of capturing the nonlinear dynamics of a structure with both types of nonlinearity.
Two categories of reduction method can be considered: linear methods and nonlinear methods, in the sense that in nonlinear reduction methods, new variables are nonlinearly related to the initial ones, while the linear methods are based on a linear change of basis [5].Using linear methods, the Reduced Order Models (ROMs) are constructed by the Galerkin projection of the equation of motion on a low-dimensional basis.These techniques are widely used in the literature due to their ease of implementation and the fact that the linear projection preserves the general expression of the equation of motion, allowing the application of the same solution strategies as those used for the full-order model.
The linear techniques face two challenges.The first is finding the best reduction basis and the second is computing the generalized nonlinear forces, which are the projection of the nonlinear forces onto the reduction basis.This work aims to address the first problem by providing the best basis in terms of the computing cost and its ability to present the nonlinear dynamics accurately.The approaches to address the second challenge in geometrically nonlinear structures and structures with friction contact are completely different.
Geometric nonlinearity is distributed, which means that the degrees of freedom (DoFs) of a significant portion (or potentially the entirety) of the structure are subjected to nonlinear elastic forces and nonlinearly coupled [5].As a result, the full displacement field is required to compute the elastic forces and a projection of the elastic forces has to be performed during the online stage at each iteration of the governing equations solution.This intrusive computation of nonlinear elastic forces, which involves nonlinear Finite Element formulation, is particularly time-consuming for a large system and results in a significant rise in the online computing cost.Several strategies are addressed to solve this issue [6][7][8].However, these remedies will be successful if the reduced space basis is selected properly.This study uses the intrusive method of computing reduced nonlinear elastic forces and focuses on finding the best reduction basis.
Due to the local nature of contact nonlinearity, a widely used method for their order reduction involves partitioning the structure into linear and nonlinear parts [9][10][11][12], where the nonlinear part includes the DoFs of the contact region (contact DoFs).The resulting reduced model keeps the contact DoFs in reduced coordinates.It eliminates the need to project the nonlinear contact forces during the online stage.
The Component Mode Synthesis (CMS) method [13], first developed for the substructuring of a linear system formed by multiple substructures, is widely used for model order reduction of structures with friction contact [9][10][11][12].The CMS technique involves three steps: decomposition of the structure into non-overlapping components, reducing the size of each component by projective reduction, and finally, coupling the individual components' ROMs through their interfaces.This allows the retention of the contact DoFs in the reduced model while replacing the other DoFs with generalized coordinates.In this way, the resulting ROM can correctly simulate the different contact states such as full-stick, gross slip, microslip, and liftoff.Generally, the DoFs kept in the reduced coordinate vector and the DoFs projected on the reduced basis are called master and slave DoFs, respectively.
CMS methods can be classified depending on the underlying reduction basis.Generally, two groups of modes are considered in CMS subspace: (i) component linear normal modes obtained by imposing specific boundary conditions at master DoFs; and (ii) static deflections due to applying either unit forces or unit displacements at master DoFs.The normal modes account for inherent component dynamics, while static deflections describe the interface behavior and are necessary for coupling the substructures.The Craig-Bampton method, as a fixed-interface-based CMS method [13], is the most widely used technique in which the static deformed shapes due to imposing unit displacement on a master DoF (constraint mode) and the normal modes of the body with fixed master DoFs are considered as reduction bases.Rubin, as a free-interface-based CMS method [14][15][16][17], is the other widely used technique.The static deformed shapes due to a unit force applied at a master DoF (attachment modes) and the normal modes of the body with free master DoFs comprise its reduction basis.
Component mode synthesis techniques have also been used for the substructuring of the systems with geometrically nonlinear substructures [18,19].Kuether et al. employed the Craig-Bampton method [18] and the accuracy of the assembled reduced models can be proved by comparing their nonlinear normal modes with those computed from the full-scale assembly.
To better model the physical behavior induced by geometric nonlinearity, Tiso and his colleagues extended CMS methods with the addition of static modal derivatives (SMDs) to the reduction basis.The formulation was first introduced in [20,21], then the extended Craig-Bampton method was more thoroughly studied in [22], in which SMDs of the constraint modes and fixed interface normal modes for each component were considered and modified as fixed-interface modes for the convenience of interface compatibility during the assembly procedure.The extension of Rubin was also studied in [23] and developed for a 3D flexible multibody system in [24], in which a nonlinear coupling between attachment modes and free-interface modes was considered.These studies concluded that extended CMS techniques can accurately predict the response, while the Rubin substructuring basis enhanced with SMDs provides a better representation of the geometric nonlinearities in the interfaces, since it comprises normal modes and SMDs featuring free motion of the interface.
Some researchers have also considered model order reduction techniques of structures with both contact and geometric nonlinearities [1,4,25,26].In [1,4,25], Craig-Bampton is employed to keep the contact DoFs in generalized coordinates.Balmaseda et al. use only the linear modes [1] while Delhez et al. augmented reduction basis with SMDs [4,25].The comparison of linear Craig-Bampton and augmented Craig-Bapton in [4] shows slightly better results of the latter while it requires a larger reduction basis to get accurate results.
To the authors' knowledge, Rubin ROM enhanced by SMDs has not yet been studied for structures with both types of nonlinearity.Given its superior performance in the substructuring of geometrically nonlinear structures, especially in capturing the effects of geometric nonlinearity at the interface [24], this model has been explored and compared to the previously used Craig-Bampton ROM augmented with SMDs.Since the number of SMDs grows quadratically with the number of modes used, an existing selection strategy [27] has been modified here to provide a more effective reduction basis.
A cantilever beam model with friction contact at its tip is considered as a case study.In this model, the accurate prediction of the coupling between transverse and axial displacements occuring in the geometrically nonlinear regime plays a critical role in correctly predicting the contact state.The Enhanced Rubin ROM is assessed to determine its ability to simulate the full-order beam behavior when both nonlinearities are active.
The formulation is outlined in Section 2, where the system of reduced-order equations required to handle a geometrically nonlinear structure with friction contact is described, along with the methods for computing the reduction basis and the projection matrix.In Section 3, the dynamics of the considered model are described and the accuracy of the aforementioned reduction methods is illustrated in different contact states.Finally, the conclusion is provided in Section 4.

Methodology
Finite Element discretization of a structure with linear elastic material undergoing large displacement/rotation and frictional contact leads to a second-order ordinary differential equation (see, e.g., Refs.[22,[28][29][30][31] for more details): in which M ∈ R N×N and C ∈ R N×N represent the mass and viscous damping matrices of the structure and u(t) ∈ R N is the vector of nodal displacement.The force vectors are the elastic internal forces f int : R N → R N , the linear excitation f e ∈ R N , and the nonlinear contact forces f c : R N → R N .The vector of internal elastic force f int accounts for both the linear and the nonlinear components.In the case of a geometrically linear structure, it is a linear function of the nodal displacement vector f int (t) = Ku(t), where K is the stiffness matrix.This study considers the mass proportional damping as C = αM.The solution of Equation ( 1) for a full-order model provides accurate results, but its computing cost is significant if the size of the system N is large.A model order reduction technique is employed to mitigate this cost by reducing the equation size.The linear projection of Equation ( 1) on a reduced space is performed by introducing the coordinate transformation: where R ∈ R N×n , the so-called projection matrix is constructed by a set of orthogonal vectors that represent the reduced space basis.Time-dependent vector q ∈ R n (n ≪ N) consists of the reduced order coordinates.Upon substitution of Equation ( 2) in Equation ( 1) and premultiplying R T in the resulting equation, the following reduced-order equation is obtained using Galerkin projection, which requires the residual to be orthogonal to the subspace span: M q(t) + C q(t) + fint (q(t)) = fe (t) + fc ( q, q, t) in which the matrices M ∈ R n×n and C ∈ R n×n represent the reduced mass and damping matrices of the structure, respectively, where Ã = R T AR.The generalized force vectors, f, are also computed as f = R T f.The linear structural matrices (M, C) and the reduced contact force f c are projected before the solution of the equations in the offline stage when the reduction basis is also computed.On the contrary, the projection of the elastic forces f int has to be performed in the online stage at each iteration of the solution of Equation ( 3) based on the predicted full displacement field.

Projection Matrix
The CMS methods start by splitting the nodal displacement vector u(t) into entries related to master DoFs u m ∈ R N m and slave DoFs u s ∈ R N s , where N m + N s = N.The master DoFs, denoted by subscript m, correspond to N m physical DoFs retained in the generalized coordinates.Other DoFs are considered slave DoFs and are denoted by subscript s.Here the master DoFs include the DoFs that the contact forces are applied to.The equation of motion Equation ( 1) and coordinate transformation matrix Equation ( 2) are then partitioned accordingly.Thus, Equation ( 2) is written as follows: where I ∈ R N m ×N m is an identity matrix.η ∈ R n g is the time-dependent vector of the amplitude of normal modes (n = N m + n g ).Accordingly, the projection matrices of Rubin (R Rubin ) and Craig-Bampton (R CB ) model order reduction techniques are computed based on the considered normal modes and static deflections ( [13][14][15]32]): Matrices Ψ r ∈ R N×N m and Φ f ∈ R N×n g , in Rubin projection matrix R Rubin , are, respectively, the matrices of the residual flexibility attachment modes and the selected normal modes with free boundary condition at master DoFs.Here, for the sake of simplicity, it is assumed that the structures are constrained and have no rigid body motion, so the rigid body modes are not considered in the Rubin basis.Ψ c ∈ R N s ×N m and Φ c ∈ R N s ×n g matrices in Craig-Bampton projection matrix R CB include constraint modes and the selected normal modes with fixed boundary condition at master DoFs, respectively.
To compute these static deflections and normal modes, first, the tangential stiffness matrix of the structure in static equilibrium (t = 0) is computed: here, f int (u(0)) = 0 and K comprise the undeformed linear stiffness matrix.Then, to obtain the free interface modes of the Rubin model, the nonlinear normal modes are computed as Equation (8).The attachment modes G m are the static displacements of the structure due to unit force applied at master DoFs.The residual flexibility attachment modes Ψ r , are computed by enforcing the spectral orthogonality between the attachment modes and the normal modes as Equation ( 9): where ω p is the natural frequency of p th free interface normal mode ϕ f p .The fixed interface normal modes and constraint modes of the Craig-Bampton method are computed from the partitions of the tangential stiffness matrix as Equations ( 10) and (11).The constraint modes Ψ c are the static displacements of the structure due to unit displacement applied at master DoFs.
In this context, the residual flexibility attachment modes and constraint modes employed in the Rubin and Craig-Bampton bases are referred to as static modes.
To account for the dependency of normal modes of a geometrically nonlinear structure on its state, the linear basis could be enhanced by adding Static Modal Derivatives (SMDs).This yields enhanced projection matrices [21,22]: where Θ f and Θ c are matrices of the selected SMDs obtained by introducing a finite displacement around the free interface and fixed interface modes, respectively.The reduction methods are characterized by the projection matrices of Equations ( 12) and ( 13) referred to as Enhanced Rubin and Enhanced Craig-Bampton methods, respectively.In the next section, the computation of SMDs is discussed.

Static Modal Derivatives
To compute SMDs, the generalized eigenvalue problem in which the nonlinear elastic force is linearized and is cast in the form f int = Ku, is differentiated with respect to the amplitude of j th mode and evaluated around the system-static equilibrium position.Ignoring the inertia terms leads to the SMDs that can also be obtained through the differentiation of the nonlinear static equation [21,22,24] as follows: where θ f ij and θ c ij are, respectively, SMDs of Rubin and Craig-Bampton bases vectors.K ′ is the derivative of the tangent stiffness matrix (K) with respect to a displacement in the direction of X j .It is prominently approximated nonintrusively via finite differences.The central finite difference yields the following formulation [33]: The SMDs do not feature orthogonality to linear modes, so a Gram-Schmidt process is employed to provide a linearly independent basis.The size of a complete set of SMDs (Θ) that can be computed from n modes is quadratic to n.However, due to their symmetry, N d = n(n + 1)/2 unique SMDs can be computed.
The number of SMDs grows quadratically with the number of linear modes.So, it is important to initially find the modes with higher participation in the structure response and calculate their corresponding SMDs.Subsequently, the most relevant SMDs are selected.To achieve this, the following steps are considered as shown in Figure 1.The suggested steps have an offline computational cost, but this procedure leads to a high accuracy through the proper selection of the SMDs and to a faster convergence.

Dominant Mode Selection
The dominant modes have a high participation factor during the system dynamics.A frequency-based approach is generally used, selecting those normal modes whose natural frequencies are close to the excitation frequency (ϕ = [ϕ (1) , . . ., ϕ (m) ]).This provides the normal modes of the linear system excited by external forces.However, the internally resonant, geometrically coupled modes, and static modes describing the local flexibility of the contact region, can contribute considerably to the behavior of a structure, particularly in strongly or even moderately nonlinear cases.Their contribution can be significant enough that the modes' derivatives are necessary in reduction bases to obtain accurate responses.The approach proposed here enables the selection of these modes.It yields a more robust reduced model that remains valid under higher levels of geometric nonlinearity.Furthermore, in terms of the convergence analysis needed to ensure the proper reduction basis, faster convergence is achieved, resulting in a shorter computing time.
To achieve this, a geometrically linear transient analysis with frictional contact (GLF analysis) along with a series of geometrically nonlinear static analyses (GNLS analyses) are conducted.The GLF analysis identifies potentially internally resonant modes excited by higher harmonic components of friction forces, while the GNLS analyses reveal geometrically coupled modes that participate significantly in the structure response.
First, the geometrically linear model with friction contact is analyzed, and the beam deflections during the steady state are maintained at different times U t = [u 1 t , u 2 t , . . ., u n t ].Then, the static force necessary to obtain a deflection equivalent to the linear dynamic response amplitude is computed as follows: when X and β are, respectively, the amplitude and the phase of the forced response of the linear system.Then, a set of nonlinear static analyses is performed at different force levels in the range [0 f] to activate different levels of geometric nonlinearity, collecting the resulting beam deflections Here, the dominant modes are defined as those with considerable participation in basis U.In other words, when U is projected onto these modes, they exhibit the highest participation factors.So, the participation of a set of modes (Υ) in displacement vectors of matrix U is computed by least square approximation U ≈ ΥQ.Then, the modes with the highest values in matrix Q will be selected as the dominant modes.This set of modes includes the free interface normal modes and residual flexibility attachment modes for the Rubin method and the fixed interface normal modes and constraint modes for the Craig-Bampton method, As shown in Figure 1, when the dominant modes Υ D are selected, the SMDs of these modes (Θ) are computed.In the next section, the selection of the most relevant SMDs is described.

Static Modal Derivative Selection
After the computation of SMDs, the next important step is the selection of the most relevant ones.The selection criterion is based on the Maximum Modal Interaction (MMI) criterion proposed by Tiso [27] for reduced-order models of geometrically nonlinear structures.This selection criterion is modified to provide higher accuracy and is adopted for structures with geometric and friction nonlinearities.The modified strategy allows the selection of the static derivatives of the static modes in case they are excited because of higher harmonics of the nonlinear forces or geometric coupling (the modes that are excited due to nonlinearity).
In its original formulation, the MMI criterion requires that a linear run under external loading over time [0, T] is first performed, and then the weighting of each SMD is computed as [20]: where η i (t) represents the time-varying amplitude of the ith mode.The SMDs with the highest value of W ij will be selected.In such a manner, only the derivatives of linearly excited modes are selected, and nonlinearity is treated as a second-order effect [20] neglecting higher orders.
Here, the matrix Q, which was previously used for the selection of the dominant modes, is employed.Matrix Q is a matrix of generalized coordinates Q = [q 1 , q 2 , . . ., q n v ] where q jT = [η j 1 , η j 2 , . . ., η j m ], n v is the number of vectors in U, and m is the total number of selected dominant modes.The weighting factor of SMDs θ ij will be computed using a similar formulation to Equation ( 19), as follows.
A modal derivative (θ ij ) whose weighting factor W ij exceeds a certain threshold is selected and added to the reduced space.Here, the weighting factors are all normalized to the largest weighting factor (W = [W 11 , W 12 , . . .]/max(W)), and a threshold of 0.0001 is considered.Unlike the original MMI approach, in this manner, the participation factors of the modes excited by nonlinear forces, such as internally resonant modes or geometrically coupled modes, are non-zero, allowing their derivatives to be selected.This results in considering the higher-order effects of nonlinearity.As a higher level of geometric nonlinearity is activated, the superior accuracy of this approach over the original MMI becomes increasingly evident.

Time Domain Analysis
The implicit Newmark method (average acceleration method) [34] with adaptive time stepping [35] is employed to obtain the time history of the transient response.In each time step, the second-order nonlinear Equation ( 3) is solved by the Newton-Raphson algorithm.This procedure is shown in Figure 2 for a geometrically nonlinear structure with friction contact.In the Newton-Raphson solver of n-th time step, as shown in Figure 2, the Newmark Beta formulation is used to obtain the first and second time-derivatives of generalized coordinate ( q(k) n+1 , q(k) n+1 ) at k-th iteration.Then, the contact forces are computed using a 2D contact model [36].To obtain the reduced nonlinear elastic forces fint (q (k) n+1 ), the physical DoFs are computed from the generalized coordinates (u = Rq), then the full-scale nonlinear elastic forces are obtained using the total Lagrangian Finite Element formulation [37].Through the projection of full-scale vector f = R T f, the reduced nonlinear elastic forces are calculated and replaced in the Newton-Raphson equation of residuals.When the residual is smaller than a defined tolerance, the Newton-Raphson solver is converged at this time step and the Newmark formulation is used to proceed in time.

Numerical Analysis
To investigate the performance of the aforementioned methods, the dynamic response of a geometrically nonlinear cantilever beam under a periodic excitation is studied.As shown in Figure 3, a friction contact is considered at its tip.The beam motion is described in plane (e x , e y ) which are the axial and transverse coordinates, respectively.Axis e x is along beam neutral axis in undeformed configuration while e y is perpendicular to it.Notably, this model can be assumed as a simplified model of blade-casing rubbing interaction in gas turbines with an assumed rigid casing and disk.The beam length is L = 2.54 m and internal structural damping is introduced as massproportional viscous damping, where α = 5.The beam material properties are Young's modulus E = 207 GPa, density ρ = 7801 kg/m 3 , and Poisson's ratio of ν = 0.28.The beam cross section is square with area A = 0.6451 m 2 .A harmonic force with the amplitude of F e = 111.29 N in a constant direction (y coordinate) is applied on the beam.A FE beam model based on a total Lagrangian nonlinear formulation and a Timoshenko kinematics [37] is considered.This model is consist of 20 elements while each node has 3 DoFs.
Besides the dynamic load applied on the contact, the contact state is defined by the static normal load N 0 and contact parameters, including contact stiffness and the friction coefficient.The effects of these parameters have been studied extensively in the literature [38].Here, the friction coefficient and contact stiffness are considered as µ = 0.5 and k t = k n = 10 N/mm.
To provide a better understanding of the dynamics of the structure under study, initially, the forced response analysis of its full order model is studied.The frequency response functions of the beam at four values of the static normal preload (N 0 = 0, 40, 60, and 250 N) are shown in Figure 4.In the case of blade-casing rubbing interaction, the normal preload N 0 is due to the centrifugal force acting on a rotating blade.The reference system is the so-called open contact (N 0 = 0 N) corresponding to the cantilever beam with no friction contact.It can be seen in Figure 4 that by increasing the value of the normal preload, initially (N 0 = 40 and 60 N) only the resonance amplitude decreases, while the resonance frequency remains constant.With the higher value of normal preload (N 0 = 250 N), the resonance frequency also increases due to the stiffening effect of the contact.
To better assess the performance of the enhanced reduction methods, two completely distinct system configurations are analyzed, corresponding to different contact states, which are representative of possible configurations of a vibrating beam potentially rubbing against a rigid casing.The name of the configuration corresponds to the contact states that dominate the nonlinear dynamics during vibration.In detail: -Slip-liftoff configuration with N 0 = 40 N and excitation frequency Freq = 3.2 Hz; -Stick-slip configuration with N 0 = 250 N and excitation frequency Freq = 8.1 Hz.

Slip-Liftoff Configuration
The dynamic behavior of the beam at the steady state when excited at the frequency of 3.2 Hz under a normal preload of 40 N is depicted in Figure 5. Figure 5a shows the beam deflections at different time instants, Figure 5b,c show the contact normal and tangential displacements, respectively, that are indeed the axial and transverse displacements of the beam tip. Figure 5b highlights that the frequency of the axial vibration is twice the frequency of excitation.In this configuration, the amplitude of vibration is 21% smaller than the maximum amplitude of the reference system (open contact).The reason for this smaller amplitude can be seen in Figure 6.The contact node state is shown in Figure 6a, where slip, stick, and liftoff are depicted by 0, 1, and −1, respectively.It shows that the contact nodes alternate liftoff and slip.As shown in Figure 6b,c, the contact tangential and normal forces are zero during liftoff.The energy dissipated by friction during slip has led to an increase in the damping ratio and a decrease in the vibration amplitude.A comparison is also made between the beam responses obtained by geometrically linear (GL) and nonlinear (GNL) models under the same preload and excitation frequency.The beam tip vibration amplitudes in transverse (U y ) and axial (U x ) directions are presented in Table 1, clearly demonstrating the inaccuracy of the GL model.As shown in Figure 7, the GL model is not capable of predicting the contact liftoff state under the current load and frequency.Consequently, a longer slip state (a higher damping) and a longer stick state (a higher stiffness) are predicted, leading to a smaller amplitude of vibration.This is due to the lack of axial-bending coupling in the GL model that underestimates the axial vibration, which plays a critical role in defining the contact state.Finally, the performance of the reduction methods is evaluated in this beam configuration.First, the SMDs selection approaches, described in Section 2.2, are investigated.Both the original MMI criterion, based on linear transient analysis, and the modified MMI criterion (MMI-NLS), based on nonlinear static analysis and geometrically linear run with friction contact, are used for the sake of comparison.
The steady-state beam response is computed by the enhanced Rubin method.Due to the importance of accurate prediction of the blade tip displacements, the error is considered to be the percentage error in the blade tip transverse vibration amplitude.This error (left axis) and the size of the reduced model (right axis) are shown in Figure 8 (left), versus the number of the normal modes used to determine the dominant modes and the SMDs.The reduced model size refers to the total number of static modes, the dominant normal modes, and the selected SMDs.Since the beam tip transverse and axial displacements (contact DoFs) are only kept as master DoFs in the reduced coordinates, there are two static modes.
The analysis is also performed with another load case (N 0 = 200 N, F e = 333.3N, and Freq = 3.2 Hz).In this case, assuming that the structure remains in the elastic region of the stress-strain curve, the contact again undergoes only slip and liftoff states and the maximum vibration amplitude is 418 mm.The geometric nonlinearity is activated to a greater extent in this case, leading to a higher contribution of geometrically coupled modes on the beam behavior.The results are shown in Figure 8 (right).The performances of Enhanced Rubin (Rubin-MD) and Enhanced Craig-Bampton (CB-MD) methods are also compared to Rubin and Craig-Bampton approaches.In particular, both SMD selection strategies are applied to the Enhanced Craig-Bampton ROM, as was previously the case for the Enhanced Rubin ROM.The results are shown in Figure 9.The 8th and 13th free interface normal modes and the 9th and 13th fixed interface normal modes are axial modes that participate in the response due to their geometric coupling with transverse modes.As seen in Figure 9 (left), adding them to the Rubin and Craig-Bampton bases considerably improved these models' accuracy.However, these reduced models with a size of 16 still exhibit a 3% error.
When the SMDs are added to the Rubin basis, the predicted response's accuracy improves significantly.Even when only the first normal mode and the two static modes (there are two contact DoFs) are considered and SMDs are computed, better accuracy is achieved than when employing large ROMs without SMDs.
Figures 8 and 9 also showcase the superior performance of the MMI-NLS approach, highlighting the relevant contribution of the selected SMDs.This is more obvious in the case of higher load (Figure 8 (right)), where the response error of Enhance Rubin ROM with three normal modes is about 0.25 % using the MMI-NLS approach and 36% using the MMI approach.The SMD of the first attachment mode (an axial mode) is the mode in which its addition to the reduction basis leads to the priority of MMI-NLS approach.This SMD is not selected by the MMI selection criterion.
To better explain the reason for the effect of the modal derivative selection on the accuracy, Table 2 shows the beam tip response error with different modes kept in the reduced space.1V refers to the first linear normal mode and 1S refers to the first static mode.In the same way, the modal derivative 1V1S is the derivative of the 1st normal mode to the first static mode.This table again shows the importance of the selection criteria.The derivative of static modes will not be selected by MMI, since these modes have no considerable participation factors during a linear run.Indeed, the SMDs of the modes that are not excited in a linear run will never be selected by MMI criterion.The resulting error cannot be ignored when geometrically coupled modes or internally resonant modes have considerable participation in the response (stronger nonlinearities).
The deformation of modal derivatives of the free interface vibration modes and residual flexibility attachment modes are shown in Figure 10, when the first free interface mode (1V) is the first bending mode.To better exhibition of axial and transverse modes, the axial displacement u x and transverse displacement u y of each node are depicted separately in the right and left plots, respectively.Figure 10 shows that the first SMD, 1V1V, is an axial mode (u y is zero in left plot).This SMD will be consistently selected via all criteria because 1V (the first transverse mode) is linearly excited by the applied harmonic force and plays a significant role in the response obtained by a linear run.Meanwhile, 2S is also a transverse mode and 1V2S is an axial mode that could be selected by the MMI criterion considering a very small threshold.On the contrary, 1S is an axial mode, which is geometrically coupled to the bending modes; 1V1S is a transverse mode (u x is zero in Figure 10 (left)) that will never be selected by the MMI approach, since the contribution factor of 1S in a linear run is zero, leading to a zero weighting factor W 1V1S = 0.This modal derivative, which is orthogonal to the first bending mode, plays a critical role in providing the best accuracy.
The deformation of modal derivatives of the fixed interface vibration modes and constraint mode is shown in Figure 11, when the first fixed interface mode (1V) is the first bending mode.The first and second constrain modes are, respectively, axial and transverse modes.They are the same as free interface modes, the static derivative of axial modes which participate due to their coupling to transverse modes; these will not be selected by MMI criterion, resulting in higher numbers of errors.

Stick-Slip Configuration
In this section, the performance of the reduction techniques is assessed in the stick-slip configuration (N 0 = 250 N and excitation frequency = 8.1 Hz).The beam deformation at different time instants during the steady state is shown in the Figure 12a, while the contact normal and tangential displacements are shown in Figure 12b,c.In this configuration, the beam amplitude of vibration decreased by about 90% compared to the reference system due to the higher normal load, leading to alternate stick and slip states.This is depicted in Figure 13, which shows the contact states and the contact forces with respect to the corresponding displacements.This contact configuration determines a higher amount of dissipated energy (see the area of the hysteresis cycle) and also an increase in the system stiffness.As a result, as illustrated in Figure 12a, not only are the resonance frequency and amplitude dependent on the contact state, and in particular on the contact stiffness, but the mode shapes of the contributing modes on structure behavior are also influenced.The beam tip vibration amplitudes obtained using the geometrical linear and nonlinear models are listed in Table 3.In this case, the linear model provides more accuracy compared to the results of the GL model in the slip-liftoff configuration, showing that the structure is in a geometrically linear regime.Although Table 3 shows that under this condition, geometric nonlinearity is not activated and the structure behavior is only influenced by the nonlinear contact forces, the accuracy of the reduction methods is compared in Figure 14. Figure 14 shows that the performance of Rubin and Craig Bampton is acceptable due to being in a geometrically linear regime.However, the enhanced Craig Bampton model has the highest accuracy with its first normal mode, two constraint modes, and five modal derivatives, which shows the effectiveness of the modal derivatives when it comes to reducing structures with friction contacts.

Conclusions
The Rubin CMS technique, enhanced with the static modal derivatives, has demonstrated its performance for the time domain analysis of a beam subjected to both local and distributed nonlinearities.The Enhanced Craig-Bampton method has also exhibited a comparable performance, particularly in a stick-slip configuration.These enhanced CMS techniques owe their performance to the addition of the static modal derivatives that allow ROMs to predict the geometrical nonlinear effects and to maintain contact DOFs that allow the accurate prediction of contact forces.The proposed strategy for the selection of modal derivatives based on the nonlinear static analysis and geometrically transient analysis with friction contact demonstrates its accuracy.The higher accuracy is due to the selection of the geometrically coupled modes and internally resonant modes as dominant modes that lead to the computation and selection of their derivatives.
Geometric nonlinear models provide more accuracy even when dealing with lowamplitude vibrations, particularly in cases when the coupling between modes could significantly impact the contact state.In the beam model with friction contact at its tip, the axial-bending coupling leads to liftoff states that cannot be predicted using a geometrically linear model.This effect can be extended to 3D models in which the coupling between the lower and higher modes is more complicated.The coupling by itself may cause small differences in structural deformation, but in the presence of friction contact, these discrepancies can result in different contact behaviors and ultimately lead to entirely different structural responses.
The nonintrusive implementation of these approaches is of particular interest, since it enables the accurate and efficient computation of the nonlinear elastic forces based on the generalized coordinates, while the distinct contact states and contact forces are anticipated accurately.

Figure 1 .
Figure 1.Selection of modes to obtain projection basis.

Figure 3 .
Figure 3. Model of a planar cantilever beam with friction contact at its end.

Figure 4 .
Figure 4. Frequency response function of the beam tip with four values of contact normal preload.

Figure 8 .
Figure 8. Response error and the size of reduction basis vs. the number of normal modes in slip-liftoff configuration (two load cases).

Figure 9 .
Figure 9. Accuracy of reduced models in response prediction to the size of reduced space in slip-liftoff configuration.

Figure 10 .
Figure 10.The transverse displacement (left) and the axial displacement (right) of the selected SMDs of free interface vibration modes and residual flexibility attachment modes.

Figure 11 .
Figure 11.The transverse displacement (left) and the axial displacement (right) of the selected SMDs of fixed interface vibration modes and constraint modes.

Figure 12 .
Figure 12.(a) Beam deformations, (b) contact normal displacement, and (c) the beam tip displacement in stick-slip configuration.

Figure 13 .
Figure 13.(a) Contact node states, (b) contact tangential, and (c) contact normal forces in stick-slip configuration.

Figure 14 .
Figure 14.Accuracy of reduced models in response prediction to the number of normal modes (N 0 = 250 N and Excitation Frequency = 8.1 Hz).

Table 2 .
Selection modes' effects on the response accuracy in slip-liftoff configuration.

Table 3 .
Vibration amplitude obtained by geometrically linear and nonlinear models with N 0 = 250 N at Freq = 8.1 Hz.