Elise Delhez

Florence Nyssen

Jean-Claude Golinval

Alain Batailly

Reduced order modeling of blades with geometric nonlinearities and contact interactions

Keywords: Reduced order model, geometrical nonlinearities, rotor/stator interaction

This article presents a methodology dedicated to the vibration analysis of turbomachine blades accounting for both geometric nonlinearities and nonlinear blade-tip/casing contacts in a numerically efficient way through the use of reduced order models. Contact is numerically handled with Lagrange multipliers and the equation of motion is integrated forward in time using an explicit central difference time integration scheme. Nonlinear internal forces caused by large displacements are evaluated using the stiffness evaluation procedure. Three reduction techniques are compared in this article, namely a nonlinear extension of the Craig-Bampton method, the proper orthogonal decomposition and a modal derivatives-based approach. These numerical methods are applied on an open industrial compressor blade model, the NASA rotor 37 blade, in order to promote reproducibility of results. The reduction methods are first applied to the blade subjected to a harmonic excitation, without contact interactions. An indicator accounting for both local and global comparison criteria is defined to more easily compare the performance of each method. When also accounting for blade-tip/casing contacts, the presented results underline that the modal derivatives-based approach is particularly well-suited for an accurate description of the blade's dynamics. This method is then employed for an in-depth analysis of the NASA rotor 37 vibration response over a wide angular speed range considering a typical contact scenario. Obtained results feature time and frequency domain signals along with stress fields. They suggest that accounting for geometric nonlinearities has a strong impact on the prediction of contact initiated interactions.

Résumé

Cet article présente une méthodologie pour l'étude des vibrations des aubes de turbomachines tenant compte des non-linéarités géométriques et des interactions de contact entre les aubes et le carter. Cette méthodologie est basée sur l'utilisation de modèles réduits afin de permettre un traitement numériquement efficace des non-linéarités. Le contact est traité numériquement à l'aide de multiplicateurs de Lagrange et les équations du mouvement sont intégrées dans le temps en utilisant un schéma d'intégration explicite par différence centrée. Les forces internes non linéaires causées par les grands déplacements de la structure sont évaluées avec la méthode dite 'stiffness evaluation procedure'. Trois méthodes de réduction sont comparées dans cet article : une variante non linéaire de la méthode de Craig-Bampton, la 'proper orthogonal decomposition' et une approche basée sur l'utilisation de dérivées modales. Ces méthodes numériques sont appliquées à un modèle ouvert d'aube de compresseur industriel (l'aube rotor 37 de la NASA) afin de permettre la reproductibilité des résultats. Les méthodes de réduction sont, dans un premier temps, appliquées à l'aube soumise à une excitation harmonique, sans interaction de contact. Un indicateur basé sur des critères de comparaison locaux et globaux est défini pour comparer les performances de chaque méthode. Dans un second temps, l'étude des interactions de contact entre l'aube et le carter montre que l'approche de réduction basée sur le concept de dérivées modales est particulièrement bien adaptée pour décrire précisément la dynamique de l'aube. Cette méthode est ensuite utilisée pour étudier la réponse vibratoire de l'aube rotor 37 de la NASA sur une large plage de vitesses de rotation, en considérant un scénario de contact typique. Les résultats obtenus dans le domaine temporel et le domaine fréquentiel ainsi que les champs de contraintes sont présentés. Ces résultats suggèrent que tenir compte des non-linéarités géométriques a une influence importante sur la prédiction des interactions de contact.

Mots-clés

Modèle réduit, non-linéarités géométriques, interaction rotor/stator 1 Introduction

In the current economic and environmental context, the main goal of aircraft engine designers is to reduce fuel consumption. From an aerodynamic point of view, one way to increase turbomachines efficiency is to reduce operating clearances between rotating blades and the surrounding casings in order to avoid air recirculation at blade tips [START_REF] Sakulkaew | Compressor efficiency variation with rotor tip gap from vanishing to large clearance[END_REF][START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF][START_REF] Nyssen | Experimental and numerical characterization of a ceramic matrix composite shroud segment under impact loading[END_REF]. As a consequence, contact interactions between rotating and stationary components are more likely to occur, even in nominal operating conditions. Today, the preliminary design of industrial turbomachine blades is mainly driven by aerodynamic considerations. There is however an increasing demand for the development of tools allowing to account for nonlinear structural considerations from the beginning of the design process in order to reduce development time.

Different numerical strategies have been developed to model and predict contact interactions [START_REF] Petrov | Generic friction models for time-domain vibration analysis of bladed disks[END_REF][START_REF] Millecamps | Snecma's Viewpoint on the Numerical and Experimental Simulation of Blade-Tip / Casing Unilateral Contacts[END_REF]. Contact can be numerically managed with penalty methods [START_REF] Sinha | Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contactimpact rub loads[END_REF], with Lagrange multipliers [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] or with the augmented Lagrangian approach [START_REF] Ma | Vibration response analysis of blade-disk dovetail structure under blade tip rubbing condition[END_REF]. In the field of turbomachinery, blade tip/casing contacts may yield several types of interactions, each of them requiring a specific numerical framework. For instance, when studying the modal interaction between the full bladed disk and the casing, the flexibility of the casing has to be accounted for in the simulations [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: Analysis for flexible casings with friction and abradable coating[END_REF]. The study of whirl motions, that can arise in accidental configurations, requires the modeling of the full bladed disk, the shaft and the casing [START_REF] Parent | Whole Engine Interaction in a Bladed Rotor-to-Stator Contact[END_REF][START_REF] Salvat | Two-dimensional modeling of unilateral contact-induced shaft precessional motions in bladed-disk/casing systems[END_REF]. When focusing on rubbing phenomena, the casing is usually supposed to be rigid but a very precise modeling of the blade/casing interface is required [START_REF] Williams | Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching method[END_REF][START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF][START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF].

In practice, numerical models of industrial bladed structures, such as finite element models, are characterized by a high number of degrees-of-freedom. The dynamic analysis of such systems therefore requires high computational capabilities. This is particularly true for nonsmooth contact simulations. Many efforts are nowadays dedicated for building reduced order models which offer a more computationally efficient framework for such numerical simulations. There exist at least two different approaches to build reduced order models. One is the so-called goal-oriented reduction approach [START_REF] Bui-Thanh | Goal-oriented, model-constrained optimization for reduction of large-scale systems[END_REF] in which the entire finite element model is directly reduced to the cost function required for an optimization. The cost function does not contain any information about the physical properties of the finite element model. This is opposite to the physical reduction approach in which the reduced order model preserves the physics of the underlying problem. This second approach has been selected in this work in order to ease the treatment of contact in the reduced space.

Among the existing physical reduction methods, a widely used one is the projection on a reduced space. It is often used in the literature because of its ease of implementation, but also because it keeps the general expression of the equation of motion so that the same integration algorithms as those used for the full order model can be used. In general, classical linear reduction methods, based for instance on the Craig-Bampton or Craig-Martinez reduction bases [START_REF] Craig | Free-interface methods of substructure coupling for dynamic analysis[END_REF][START_REF] Seshu | Substructuring and Component Mode Synthesis[END_REF][START_REF] Craig | Coupling of substructures for dynamic analyses: An overview[END_REF], can be used to reduce the system while providing accurate prediction of contact dynamics [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers : Analysis for rigid casings[END_REF][START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF]. In the case of strong contact nonlinearities, the nonlinear nature of the system may be included in the reduction basis, for instance by the use of nonlinear complex modes [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF][START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF].

The reduction methods generally used to study the contact dynamics of blades do not account for the possible large displacements and large deformations of the structure. This has not been a real issue so far as contact interactions mainly occur in compressor stages while slender blades undergoing large displacements are essentially located in fan or turbine stages. Nowadays however, all blades are designed increasingly lighter and more flexible to be in line with the tight weight constraints in aeronautics [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Grolet | Vibration Analysis of a Nonlinear System With Cyclic Symmetry[END_REF][START_REF] Picou | Mistuning analysis of a detuned bladed-disk with geometrical nonlinearities[END_REF]. Large displacements and large deformations of the blades, leading to a geometrically nonlinear behavior of the structure, can no longer be neglected [START_REF] Xiao | Interactions between blades and abradable coatings: a numerical approach considering geometrical nonlinearity[END_REF].

The model order reduction of mechanical systems undergoing large displacements started to gain the attention of researchers in the last few years. The different projective reduction methods differ by the choice of the reduction basis and by the way the nonlinear internal forces are evaluated in the reduced space. Different projection bases have been investigated in a nonlinear context: linear bases [START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF][START_REF] Palma | Nonlinear Harmonic Analysis of a Blade Model Subjected To Large Geometrical Deflection and Internal Resonance[END_REF], bases obtained from Proper Orthogonal Decomposition (POD) [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: An overview[END_REF][START_REF] Picou | Mistuning analysis of a detuned bladed-disk with geometrical nonlinearities[END_REF] and linear bases augmented with modal derivatives [START_REF] Weeger | On the use of modal derivatives for nonlinear model order reduction[END_REF][START_REF] Wu | Nonlinear model order reduction for flexible multibody dynamics: a modal derivatives approach[END_REF][START_REF] Karamooz Mahdiabadi | A non-intrusive model-order reduction of geometrically nonlinear structural dynamics using modal derivatives[END_REF]. The reduced nonlinear forces can be identified with intrusive and non-intrusive methods. In intrusive methods, reduced nonlinear forces are obtained by direct projection of the full order nonlinear forces [START_REF] Capiez-Lernout | Computational stochastic statics of an uncertain curved structure with geometrical nonlinearity in three-dimensional elasticity[END_REF][START_REF] Touzé | Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models[END_REF][START_REF] Chaturantabut | Discrete Empirical Interpolation for Nonlinear Model Reduction[END_REF]. These methods are based on the relationship between the finite element tensors and their modal counterparts. However, the expression of the full order nonlinear forces cannot always be extracted from commercial finite element software packages. Non-intrusive methods compatible with standard finite element software packages have therefore been derived. In such methods, the reduced nonlinear Delhez et al.

forces are expressed as a polynomial form of the reduced coordinates whose coefficients are identified with nonlinear static evaluations with imposed displacements (enforced displacement method or STiffness Evaluation Procedure (STEP) [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF]) or imposed forces (applied force method [START_REF] Mcewan | A combined modal/finite element analysis technique for the dynamic response of a non-linear beam to harmonic excitation[END_REF]).

To the best of the author's knowledge, none of the numerical methods mentioned above have ever been used in combination with structural contacts related nonlinearities. Thus, there is a need to assess the applicability of these methods to the analysis of the combined effect of geometric and contact nonlinearities on industrial structures.

The objective of this article is to show that the numerical strategies used for the modeling of contact interactions in turbomachines can be combined with those used for the modeling of geometric nonlinearities in order to study in a computationally efficient way the vibration response of blades subjected to contact interactions while accounting for geometric nonlinearities. In particular, this paper focuses on non-intrusive projection-based reduction methods to be compatible with commercial finite element software. Computations are made in the time domain in order to capture transient phenomena, which seem key for the rise of some critical interactions [START_REF] Millecamps | Influence of Thermal Effects During Blade-Casing Contact Experiments[END_REF]. It is also checked that the reduction methods are able to predict the steady-state dynamics. Given the complexity of the numerical simulations and the number of nonlinear degrees-of-freedom, the use of frequency methods is nowadays not really feasible for this kind of problem.

The paper is organized as follows. The specific context of the research is first recalled in section 2, together with the description of the test case analyzed in the paper. Then, three different reduction methods are selected as interesting candidates for the subsequent studies, respectively the POD, the Craig-Bampton method with POD filtering of the nonlinear forces, and a modal derivatives-based approach. These methods are described and adapted to allow the definition of contact interfaces in the reduced space in section 3. The numerical strategy adopted for contact modeling is also detailed. In section 4, the different reduction methods are applied to an industrial test case subjected to an equivalent aerodynamic loading without contact interactions. The performance and the robustness of the different reduction methods are assessed. Lastly, the reduction methods are applied to the industrial test case undergoing large displacement due to contact interactions in section 5. The most suitable reduction method is selected and used to perform detailed analyses, focusing on rubbing interactions.

Research context and test case

Over the past years, many research works related to the analysis of turbomachine blades nonlinear vibrations have been published. Researchers have focused on a variety of nonlinearities related to: the shaft/bearing interface [START_REF] Lesaffre | Contact analysis of a flexible bladed-rotor[END_REF], the blade/disk interface [START_REF] Krack | Reliability optimization of friction-damped systems using nonlinear modes[END_REF], the blade tip/casing interface [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers : Analysis for rigid casings[END_REF] or large displacements [START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF]. While obtaining relevant results for industrial mechanical systems calls for the use of realistic full 3D finite element models, the latter often prevent full disclosure of the results for confidentiality reasons. It is thus extremely arduous to compare existing numerical strategies dedicated to such analyses. Similarly to what was successfully done in the field of computational fluid dynamics [START_REF] Qi | Establishment of an open 3D steam turbine flutter test case[END_REF], researchers recently proposed to use an open industrial blade model dedicated to nonlinear structural dynamics simulations [START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF]. The establishment of this benchmark comes with a detailed description of the blade, with all necessary geometry parameters and material properties to ease reproducibility and scientific collaboration between research teams. This industrial blade model, based on NASA rotor 37 blade, is used throughout this article.

Open blade model

The benchmark proposed in [START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF] is based on the open geometry of the NASA rotor 37 compressor blade. The rotor 37 is a transonic axial flow compressor stage composed of 36 blades (see Fig. 1). It was initially designed and tested at NASA's Lewis Research Center in the late 70's [START_REF] Moore | Performance of Single-Stage Axial-Flow Transonic Compressor With Rotor and Stator Aspect Ratios of 1.19 and 1.26[END_REF][START_REF] Reid | Design and Overall Performance of Four Highly-Loaded, High Speed Inlet Stages for an Advanced High-Pressure Ratio Core Compressor[END_REF]. Since then, NASA rotor 37 has been often used as a benchmark to validate numerical methods in the field of aerodynamics [START_REF] Denton | Lessons from Rotor 37[END_REF].

Rotor 37 is made of 200-grade maraging steel, a nickel alloy [START_REF] Reid | Design and Overall Performance of Four Highly-Loaded, High Speed Inlet Stages for an Advanced High-Pressure Ratio Core Compressor[END_REF], with Young's modulus E = 180 GPa, Poisson's ratio ν = 0.3 and density ρ = 8, 000 kg/m 3 . The 3D geometry of the blade is described with a finite element model composed of quadratic pentahedral elements. The discretized mesh of the blade is represented in Fig. 2. It counts 20,657 nodes, i.e. 61,971 degrees-of-freedom. The compressor blade is clamped at its root. The first three eigenfrequencies of the blades are given by ω 1 = 5, 691.5 rad/s, ω 2 = 17, 062.9 rad/s and ω 3 = 20, 564.2 rad/s. The Delhez et al. Single blade of NASA rotor 37 with mesh used for the finite element simulations and clamped surface ( ). The seven colored nodes are the blade tip nodes used further for the definition of the excitation and for contact simulations. Among those nodes, the node T E is located at the trailing edge (•) and the node LE is located at the leading edge (•).

Geometric nonlinearities influence

As previous studies have shown that single blade numerical models led to results in good agreement with experimental measurements when focusing on rubbing interactions, the dynamics of the disk and the surrounding blades is neglected [START_REF] Batailly | Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF]. The spatially discretized equation of motion governing the physical displacement u of the blade around its equilibrium state writes

Mü + C u + Ku + g nl (u) = f e (t) + f c (u, u), ( 1
)
where M is the mass matrix, C is the viscous damping matrix, K is the linear stiffness matrix, g nl (u) is the nonlinear internal forces vector, f e (t) is the external forces vector and f c (u, u) accounts for possible contact interactions. The time dependence of u is omitted for clarity. Regarding the modeling of damping, the assumption of proportional damping is made. As suggested in [START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF], the damping matrix C is computed in the modal domain with a modal damping coefficient ξ = 5 • 10 -4 .

Delhez et al. The reduction method used in [START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF] does not allow to account for geometric nonlinearities. However, numerical simulations performed on the finite element model of the rotor 37 blade show that geometric nonlinearities may affect significantly the response of the structure, even for small amplitudes of vibration. To illustrate this, the blade is excited at the leading edge (node LE in Fig. 2) in the e y direction by a harmonic excitation of amplitude A e = 300 N and pulsation ω e = 2, 000 rad/s. Contact interactions are not considered, i.e. f c (u, u) = 0 in Eq. [START_REF] Sakulkaew | Compressor efficiency variation with rotor tip gap from vanishing to large clearance[END_REF]. The responses of the blade at node LE in the e x direction to this excitation are compared in Fig. 4 when geometric nonlinearities are neglected (i.e. when Eq. ( 1) is linearized by taking g nl (u) = 0) and when they are taken into account. Significant differences are observed between both signals. In particular, the nonlinear response has a much richer frequency content than the linear one. This is patent in Fig. 5 where the responses are compared in the frequency domain. For the linear solution, energy is only localized at the frequency of the excitation and at the frequencies of the excited mode shapes. For the nonlinear solution, harmonics and subharmonics of the excitation pulsation and eigenfrequencies are also excited. This example shows that geometric nonlinearities may have a strong effect on the dynamics of the blade, even for displacement amplitudes smaller than 1 mm that could for instance be induced by contact interactions with the casing. It is therefore of considerable importance to be able to account for the geometric nonlinearities in the reduced order models used to predict contact interactions. This paper aims at showing how reduced order modeling techniques suitable for contact simulations can be enriched to account for geometric nonlinearities. 3 Reduced order modeling

Projection-based model order reduction

The dynamics of mechanical structures is usually studied with high fidelity finite element models. The spatially discretized equation of motion governing the physical displacement u of the structure around its equilibrium state is given by Eq. ( 1). In projection-based model order reduction techniques, the displacement field u is written as a linear combination of basis functions

Φ u = Φ q. (2)
The physical displacement u ∈ R NDOF is therefore linearly reduced to a set of generalized coordinates q ∈ R r , where r N DOF (in practice r N DOF ). Ideally, the reduced order basis Φ must be chosen such that it spans a subspace to which the solution of the full order model belongs. Injecting Eq. (2) in Eq. ( 1), the projection of the equation of motion in the reduced space writes Mq + C q + Kq + gnl (q) = fe (t) + fc (q, q), (

where the projected matrices take the expressions

M = Φ T MΦ, C = Φ T CΦ, K = Φ T KΦ (4)
and the projected vectors

gnl (q) = Φ T g nl (Φq), fe (t) = Φ T f e (t), fc (q, q) = Φ T f c (Φq, Φ q). ( 5
)
As mentioned in the introduction, this way of reducing a mechanical system advantageously preserves the general form of the equation of motion. The initial conditions on the generalized coordinates related to the reduced model are obtained from the initial conditions on the degrees-of-freedom of the full order model by a least-squares approximation, i.e.

q(0) = Φ T Φ -1 Φ T u(0) and q(0) = Φ T Φ -1 Φ T u(0). ( 6
)
The projection-based model order reduction methods differentiate between an expensive offline stage where the projection basis is computed and the reduced order model is built and an efficient online stage where the problem is solved. Figure 6 schematically represents the distinction between offline and online stages. The offline stage takes the full order finite element model as an input. In this stage, the reduction basis is built, the full order matrices are projected onto the reduction basis and the expression of the nonlinear internal forces is identified with an Delhez et al.

adequate method. In the online stage, numerical simulations are performed. In this stage, the initial conditions and the excitation are first projected onto the reduced space. The equation of motion in the reduced space is then numerically integrated in time before the solution may be projected onto the full order space. The offline stage is performed only once for a given reduction basis, while the online stage can be performed several times to study different excitation conditions. When building nonlinear reduced order models with projective reduction methods, different challenges must be addressed: (1) finding an adequate reduction basis embedding the main nonlinear characteristics of the structure and (2) evaluating the internal nonlinear forces and the contact forces in the reduced space. The next three paragraphs focus on these two challenges.

Reduction basis

Three different reduction bases have been selected as good candidates to reduce mechanical systems featuring geometric nonlinearities and contact nonlinearities:

• reduction bases obtained by proper orthogonal decomposition, which are widely used when studying the dynamics of structure [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: An overview[END_REF]; • the Craig-Bampton reduction basis, which is largely used in a linear context and has also proved to be efficient when considering contact nonlinearities [START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF]; • the concept of modal derivatives, which begins to gain the attention of researchers for the modeling of geometric nonlinearities [START_REF] Weeger | On the use of modal derivatives for nonlinear model order reduction[END_REF]. In order to keep some physical degrees-of-freedom in the reduced model, and therefore to allow the treatment of contact constraints directly in the reduced space, the different reduction bases are -when required -adapted to retain the corresponding physical degrees-of-freedom in the reduced coordinates. In particular, as detailed further, the reduction method by proper orthogonal decomposition requires some adjustments and a variant (based on the Craig-Bampton modes) of the classical modal derivative approach is selected. The physical displacement vector u is therefore partitioned into its inner and boundary degrees-of-freedom counterparts, u i and u b as

u = u i u b . (7)
The boundary degrees-of-freedom should be included in the reduced space, i.e. the degrees-of-freedom u b are part of the generalized coordinates q. The structural matrices M, C and K and the force vectors g nl , f e and f c are partitioned accordingly with respect to inner (subscript i) and boundary (subscript b) degrees-of-freedom

Y = Y ii Y ib Y bi Y bb and y = y i y b , (8)
Delhez et al.

where Y stands for the mass, damping or stiffness matrix and y stands for the nonlinear, external or contact force vector. The vector of generalized coordinates is partitioned as

q = u b η , (9)
where η stands for the modal generalized coordinates that do not have any physical meaning.

Proper Orthogonal Decomposition

The Proper Orthogonal Decomposition (POD), also known as the Karhunen-Loève decomposition, was proposed independently by several authors in the 40's [START_REF] Karhunen | Uber Lineare Methoden in der Wahrscheinlichkeitsrechnung [About linear methods in probability[END_REF][START_REF] Loeve | Fonctions aléatoires du second ordre[END_REF]. The model reduction by POD is a data-driven method which identifies the optimal subspace from a set of displacement snapshots obtained from simulations of the full order model. If the m displacement snapshots are denoted by u (t i ) (i = 1, . . . , m) and collected in a matrix

X = u (t 1 ) u (t 2 ) • • • u (t m ) , (10)
the principle of the POD reduction method consists in finding the subspace Φ such that the projection of the data matrix X onto Φ has maximal spatial variance. The POD basis can be obtained by identifying the dominant left singular vectors of the Singular Value Decomposition (SVD) of the data matrix X corresponding to the highest singular values.

In its general formulation, the POD method does not allow to retain physical degrees-of-freedom in the generalized coordinates q. In order to ease the implementation of contact interfaces, it is proposed to adapt this reduction method by defining the projection basis Φ as a combination of r b constraint modes and r POD POD modes.

The constraint modes, or Guyan modes, correspond to the linear static deformation of the structure to unit displacements at the boundary degrees-of-freedom. They can be written as

Ψ = Ψ i I , (11)
where Ψ i is the solution of

KΨ = K ii K ib K bi K bb Ψ i I = 0 R , (12)
which gives

Ψ i = -K -1 ii K ib . (13)
The POD modes are computed in the same way as described previously. According to the degrees-of-freedom partitioning, they write

Φ s = Φ s,i Φ s,b . (14)
The interface component of the POD modes is not null since the interface is not fixed. Since the interface motion is already fully represented by the constraint modes, the POD vectors are modified to fixed-interface modes by subtracting the component already present in the constraint modes

Φ POD 0 = Φ s,i Φ s,b - Ψ i I Φ s,b . (15)
The reduction basis therefore writes

Φ = Ψ i Φ POD I 0 . (16)
Delhez et al.

Craig-Bampton basis

The Craig-Bampton reduction method has been widely applied in linear contexts [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF], but also in the presence of contact nonlinearities [START_REF] Piollet | Blade/casing rubbing interactions in aircraft engines: Numerical benchmark and design guidelines based on NASA rotor 37[END_REF]. In this method, the projection basis Φ is composed of r c fixed interface linear normal modes Θ and r b constraint modes Ψ. The fixed interface linear normal modes correspond to the modes of the structure clamped at the boundary interface degrees-of-freedom. They can be written as

Θ = Θ i,rc 0 , (17)
where Θ i,rc is formed by truncating the solution Θ i of the eigenvalue problem

K ii Θ i = M ii Θ i ω 2 (18)
to the first r c modes. The constraint modes are defined by Eq. ( 11). The reduction basis therefore writes

Φ = Ψ i Θ i,rc I 0 . ( 19
)
Different studies have shown that, without correction, this linear reduction basis is not suitable for the modeling of geometric nonlinearities [START_REF] Grolet | Vibration Analysis of a Nonlinear System With Cyclic Symmetry[END_REF][START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF]. Indeed, even if some high-frequency linear modes could capture the geometrically nonlinear behavior of the structure, these modes are truncated in the reduction procedure. Including these modes in the reduction basis would lead to reduction bases of large size and the benefit of the reduction procedure would be lost [START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF][START_REF] Vizzaccaro | Nonintrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using threedimensional finite elements[END_REF]. However, as shown in the remainder, it is possible to correct the evaluation of nonlinear internal forces based, for instance, on a POD filtering approach. Such correction allows to build reduced order models accurate and robust in the nonlinear regime. Another approach is to enrich the linear bases with modes allowing to intrinsically account for the nonlinear behavior of the structures such as dual modes [START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF] or modal derivatives, as described in the next section.

Craig-Bampton basis enhanced with modal derivatives

The inherent dependence of a nonlinear mechanical system's modal properties on its total energy level implies that linear reduction methods are ill-suited in a nonlinear framework. Enriching linear bases with modal derivatives was first introduced as a way of intrinsically accounting for the nonlinear behavior of the structure in linear bases composed of linear normal modes [START_REF] Tiso | Optimal second order reduction basis selection for nonlinear transient analysis[END_REF]. The concept has then been recently extended to other linear reduction bases, and in particular to the classical Craig-Bampton basis [START_REF] Wu | A modal derivatives enhanced Craig-Bampton method for geometrically nonlinear structural dynamics[END_REF].

Let's denote by Φ(u) the Craig-Bampton modes computed at a certain level of deformation u. In practice, these modes are computed in the same way as the modes Φ in Eq. ( 19) by using the tangent stiffness matrix

K T = K + ∂g nl ∂u (20)
instead of the stiffness matrix K in Eq. ( 13) and Eq. ( 18). Modal derivatives write

ψ ij = ∂ φi ∂q j q=0 i, j = 1, • • • , r, (21)
where φi stands for the i th mode of the reduction basis Φ and q j is the j th generalized coordinate, see Eq. ( 2). Modal derivatives describe the variation of the linear modes φ i when the structure is deformed in the shape of another mode. In practice, modal derivatives are computed by solving the linear systems

Kψ ij = - ∂K T ∂q i q=0 φ j . (22)
It can be shown that modal derivatives are symmetric [START_REF] Weeger | On the use of modal derivatives for nonlinear model order reduction[END_REF], i.e.

ψ ij = ψ ji , (23)
so that a basis consisting in r linear modes can be enhanced with r(r + 1)/2 modal derivatives.

The size of the reduction basis therefore increases quadratically with the number of linear modes. However, it is expected that only a few modal derivatives are necessary to capture the nonlinear behavior of the structure [START_REF] Tiso | Effective modal derivatives based reduction method for geometrically nonlinear structures[END_REF]. In this article, a selection criterion is therefore defined in order to select the best modal derivatives. In the selection procedure, a weight is associated to all modal derivatives according to their contribution and the modal derivatives characterized by the highest weights are selected. The weight is defined by the modal interaction between the modes during a linear simulation

W ij = T 0 |q i (t)q j (t)|dt, ( 24
)
where W ij represents the weight associated to the modal derivative ψ ij and q i (t) is the time varying amplitude of the i th linear mode obtained in the linear response to a given applied external loading. This selection criterion based on the modal interaction between the modes has proven to be efficient to select modal derivatives of linear normal modes [START_REF] Tiso | Optimal second order reduction basis selection for nonlinear transient analysis[END_REF]. It is here proposed to apply this selection criterion for the selection of modal derivatives of Craig-Bampton modes. While the computation of the modal derivatives is a simulation-free procedure, i.e. a procedure that does not require the time integration of the full order equation of motion (contrary to the POD method for instance), the selection of the modal derivatives is based on the numerical time integration of Eq. ( 1) where the nonlinear terms are neglected (i.e. g nl (u) = 0 and f c (u, u) = 0).

Computation of the reduced nonlinear internal forces

The direct computation of the reduced nonlinear forces gnl (q) with Eq. ( 5) is not feasible for two main reasons. First, the computational complexity of the nonlinear term evaluation depends on the total number of degrees-of-freedom. For industrial structures, the cost of evaluating the nonlinear term increases rapidly to a point that the computational gain through the use of reduced order model is lost. Then, this way of proceeding is intrusive because it requires to have access to the full order internal forces, and therefore to the internal models of the commercial finite element software. In this paper, the reduced nonlinear term is evaluated with the STiffness Evaluation Procedure (STEP). This non-intrusive method is compatible with commercial finite element software packages and has proven to be efficient for the modeling of dynamical mechanical systems subjected to large displacements [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF][START_REF] Capiez-Lernout | Geometric nonlinear dynamic analysis of uncertain structures with cyclic symmetry. Application to a mistuned industrial bladed disk[END_REF][START_REF] Givois | On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models[END_REF]. In the case where the reduction basis does not intrinsically account for the nonlinear behavior of the structure, the STEP provides poor results but can be easily adapted to build reduced order models model accurate in the nonlinear regime [START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF].

STiffness Evaluation Procedure (STEP)

The STEP has been first introduced in [START_REF] Muravyov | Determination of nonlinear stiffness with application to random vibration of geometrically nonlinear structures[END_REF]. It consists in imposing static displacement fields and determining the corresponding internal nonlinear forces. In the STEP, the nonlinear reduced term is assumed to be a polynomial of degree 3 in terms of the generalized coordinates q 1 , q 2 , . . . , q r , components of the vector of generalized coordinates q.

The degree of this polynomial reflects the fact that the nonlinearities in the nonlinear internal forces g nl (u) are cubic under the assumption of a linearly elastic material. A polynomial of third order could actually be an exact representation of the nonlinear internal forces for the full order finite element model. However, by being applied to the reduced order system, this polynomial becomes an approximation. The polynomial formulation reads

gnl (q) = Ãq ⊗2 + Bq ⊗3 , ( 25
)
where à is a tensor in R r 2 r and B is a tensor in R r 3 r . In the following, the single elements of the tensors will be denoted by

Ãm ij = Ã([i, j] , m) and Bm ijk = B([i, j, k] , m). (26)
Single-dimension-tensors extracted from the tensors à and B are denoted by

Ãij = Ã([i, j] , :) and Bijk = B([i, j, k] , :). ( 27
)
The components of the (r × 1) vector gnl (q) are approximated by

gm nl (q) = r i=1 r j=i Ãm ij q i q j + r i=1 r j=i r k=j Bm ijk q i q j q k , (28)
with m = 1, 2, . . . , r. The limited ranges of the indices can be exploited to reduce computational and storage requirements.

The coefficients Ãm ij and Bm ijk appearing in Eq. ( 28) are determined by nonlinear forces static evaluations performed on the full order finite element model. First, the coefficients with equal indices Ãm ii and Bm iii (i = 1, . . . , r) are determined by imposing displacements of the form

u = +q i φ i and u = -q i φ i (i = 1, . . . , r), (29)
where φ i is the i th mode of the reduction basis Φ, and computing the associated nonlinear forces

F NL a = g nl (u = +q i φ i ) and F NL b = g nl (u = -q i φ i ) . ( 30
)
The absolute values of the modal coordinates q i imposed must be chosen sufficiently large to activate the nonlinear behavior of the structure but sufficiently small to remain in the convergence limit of the finite element software. The Galerkin projections of these forces provide the reduced nonlinear forces

FNL a = Φ T F NL a = Φ T g nl (u = +q i φ i ) = gnl q = [0 • • • + q i • • • 0] T (31)
and

FNL b = Φ T F NL b = Φ T g nl (u = -q i φ i ) = gnl q = [0 • • • -q i • • • 0] T . ( 32
)
Combining Equations ( 28), ( 31) and (32) yields the system of equations

FNL a = Ãii q 2 i + Biii q 3 i FNL b = Ãii q 2 i -Biii q 3 i ( 33
)
where q i , FNL a and FNL b are known. The vectors Ãii and Biii can therefore be identified

             Ãii = FNL a + FNL b 2q 2 i Biii = FNL a -FNL b 2q 3 i ( 34
)
The coefficients with two unequal indices Ãm ij , Bm iij and Bm ijj (i < j) can be identified in the same way by imposing displacements of the form u = +q i φ i + q j φ j , u = -q i φ i -q j φ j and u = +q i φ i -q j φ j (i, j = 1, . . . , r ; i < j).

(

) 35
The coefficients with three unequal indices Bm ijk (i < j < k) can be identified by imposing displacements of the form

u = +q i φ i + q j φ j + q k φ k (i, j, k = 1, . . . , r ; i < j < k). (36)
The total number N static of required nonlinear forces static evaluations is given by

N static = r 3 6 + r 2 + 5r 6 , ( 37
)
where r is the number of modes retained in the reduction basis. This number can be viewed as a measure of the fixed cost of the reduced order analysis as the model reduction must be performed regardless of the simulated response time to be computed.

POD-based correction

In the case where the projection basis is built from linear computations only (such as the Craig-Bampton basis introduced in section 3.2.2, or bases composed of linear vibration modes), the STEP provides poor results [START_REF] Grolet | Vibration Analysis of a Nonlinear System With Cyclic Symmetry[END_REF][START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF]. More specifically, the projection of nonlinear forces onto a linear basis to obtain the generalized forces in the reduced order model may create spurious artifacts in the time response, inducing low amplitude displacement with high harmonics, as the linear modes that could represent the nonlinear behavior of the structure are high-frequency modes that are generally truncated in the reduction procedure [START_REF] Mignolet | A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures[END_REF]. To avoid such behavior, it has been suggested to carry out a filtering of the full order nonlinear forces [START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF]. In practice, the nonlinear forces are collected from a given number of snapshots that represent a set of characteristic displacements in the response u 1 , ..., u m ,

X = g nl (u 1 ) g nl (u 2 ) • • • g nl (u m ) . (38)
The nonlinear basis Φ f used for the filtering is obtained by implementing the SVD of matrix X. The nonlinear basis consists in the truncation to r f modes in the resulting left singular vectors of the SVD basis. The nonlinear forces in the full order model are then filtered with this basis Φ f as

g f nl (u) = Φ f Φ T f g nl (u) (39)
and the resulting filtered forces g f nl (u) are used as an input to the STEP method instead of the exact nonlinear forces g nl (u) .

Contact algorithm

As model order reduction by projection preserves the general form of the equation of motion and the different adapted reduction techniques introduced in section 3.2 allow to retain the physical contact interface in the reduced space, contact algorithms used for full order finite element models can be directly used for reduced models.

The numerical strategy adopted here relies on the explicit central difference time integration scheme with the handling of contact using Lagrange multipliers [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF]. Such simulations in the time domain allow to capture the transient dynamics of the structure due to the occurrence of contact events. Friction is managed with the Coulomb law, under the assumption of sliding friction as the tangential speed to rotation is significantly larger than the vibration speed of the structure.

By noting q n , the numerical approximation of the exact value q(t n ), the classical explicit central finite difference scheme used in this work yields qn = q n+1 -2q n + q n-1 h 2 and qn =

q n+1 -q n-1 2h , ( 40
)
where h is the time step of integration. At each time step t n+1 = t n + h, the displacement field is first predicted without accounting for possible contact interaction as

q n+1 = 2 M + h C -1 2h 2 fe (t n ) -gnl (q n ) (41) + 4 M -2h 2 K q n + h C -2 M q n-1 .
Delhez et al.

Then, the gap between the impacting nodes and the target surface is computed. If a penetration is detected, the contact forces fc,n+1 are computed in the normal and tangential directions through the use of Lagrange multipliers using the methodology described in [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] in order to cancel the penetrations of the impacting nodes. The displacement field is finally corrected to avoid penetration by

q n+1 = q n+1 + M h 2 + C 2h -1 fc,n+1 . (42)

Analysis with geometric nonlinearities alone

As a first step, the three different reduction methods introduced in the previous section are used to predict the response of the rotor 37 blade to a harmonic excitation, representative of a typical aerodynamic excitation, without contact interaction.

Construction of the reduced order models

The reduced order models are built in the offline stage of the numerical procedure. In order to ease the comparison, the physical interface kept in the reduced space only counts one node. The r b = 3 degrees-of-freedom corresponding to the node LE at the blade tip trailing edge are kept in the reduced space. For the three reduction methods, the nonlinear internal forces are evaluated with the STEP (see section 3.3.1). The values of the generalized coordinates q i , q j and q k appearing in Eq. ( 29), Eq. ( 35) and Eq. ( 36) are defined so that the blade reaches displacement amplitudes of the order of 1 cm in the nonlinear static analyses performed in the STEP. The convergence of some elements of the tensors à and B identified in the procedure with the values of q i , q j and q k has been studied in order to check that the geometrically nonlinear regime is reached. The characteristics of the reduction bases are summarized in Tab. 1. In order to study how richer reduction bases can enhance the accuracy of the predictions, different sizes r of the reduction bases are considered for each method. POD. The first reduction basis is built by combining (r -3) POD modes and 3 constraint modes (as described in section 3.2.1). The displacement snapshots used to build the basis are obtained by computing the response of the blade to a training excitation with the full order model. The blade is here excited at node LE in the e y direction with a harmonic excitation during one period. The harmonic excitation is characterized by an amplitude of 300 N and a pulsation of 2,000 rad/s. The excitation frequency allows to mainly excite the first bending mode of the blade. The corresponding nonlinear response at node LE is shown in Fig. 4.

The comparison with the linear response confirms that the high amplitude of the excitation allows to excite nonlinearities in the structure. Craig-Bampton with POD filtering. The Craig-Bampton basis is built with 3 constraint modes and (r -3) fixed interface modes, as detailed in section 3.2.2. Without correction, this reduced model does not allow to capture the nonlinear nature of the system. This is illustrated in Fig. 7, which compares the reference nonlinear solution of the structure with the solution obtained with the usual linear Craig-Bampton reduction method (r = 15) where the nonlinear internal forces are evaluated with the STEP. Very poor results are obtained. As proposed in [START_REF] Balmaseda | Reduced order models for nonlinear dynamic analysis with application to a fan blade[END_REF] and explained in section 3.3.2, results can be enhanced by filtering the nonlinear internal forces through a POD correction before applying the STEP. The training excitation used for the POD filtering of the nonlinear internal forces is here defined as a harmonic excitation of amplitude A t = 300 N and a pulsation of ω t = 2, 000 rad/s applied at node LE in the e y direction.

Modal derivative. The modal derivative basis is built as described in section 3.2.3, by combining 3 constraint modes, 4 fixed interface modes and (r -7) modal derivatives selected according to the criteria defined by Eq. [START_REF] Xiao | Interactions between blades and abradable coatings: a numerical approach considering geometrical nonlinearity[END_REF]. Figure 8 shows the weights W ij associated to the different modal derivatives for the training excitation already used to build the POD reduction basis and the Craig-Bampton reduction basis. The weights have been raised to the power of 0.25 to better show the relative contributions of the modes. This figure also highlights the 10 modal derivatives with the highest weights (blue bars in the figure). The dominating modal derivatives correspond to the modal derivatives of the static modes with respect to the static modes themselves. Then, the modal derivatives of the static modes with respect to the first internal modes have to be selected. It should be noted that it is not guaranteed that the modal derivatives are orthogonal and linearly independent. The modal derivatives are therefore orthonormalized with a modified Gram-Schmidt procedure which also allows to remove the linearly dependent vectors from the basis. This transformation allows to get better results and to accelerate the convergence of the nonlinear time integration scheme.

Performance indicator

In order to ease the comparison of the results obtained with each method, a performance indicator combining different quantities of interest is defined. Before being used to assess the accuracy and the versatility of the reduction methods in sections 4.3.1 and 4.3.3, this performance indicator is used to check the relevance of the modal derivatives selection criterion at the end of the section.

Definition

The performance indicator I is defined is a weighted sum of N c different criteria

I = Nc i=1 w i e i . (43)
Each criterion e i varies between 0 (worst case) and 1 (perfect agreement between the solutions obtained with the full and reduced order models). The sum of the weighting factors w i equals 1. The relative magnitudes of the weighting factors reflect the relative influence of the different criteria in the performance indicator definition. The performance indicator therefore also ranges from 0 to 1, the higher the performance indicator, the better the reduction method. Four criteria are here selected in order to account for the global and local differences between the time solutions obtained with the full order finite element model u FOM (t) and the reduced order model u ROM (t) = Φ q(t).

• The first criterion is defined as the covariance criterion at maximal force amplitude

e 1 = x T 1 x 2 2 x T 1 x 1 x T 2 x 2 , (44)
where x 1 = u FOM (t max ), x 2 = u ROM (t max ) and t max corresponds to the time when the excitation reaches its first maximal value. • The second criterion is computed from the relative error on the energy (defined as the sum of the kinetic energy and the linear energy of deformation) and writes

e 2 = max   0, 1 - t [∆E(t)] 2 t [E FOM (t)] 2   , (45)
where

E(t) = 1 2 u(t) T M u(t) + 1 2 u(t) T Ku(t) (46)
and

∆E(t) = E FOM (t) -E ROM (t). (47)
• The third criterion is based on the relative oscillations deviation at blade tip and is defined as

e 3 = max   0, 1 - 1 3N tip i∈indextip o FOM i -o ROM i o FOM i 2   , (48)
where N tip is the number of nodes located at blade tip, index tip stands for the indices corresponding to the degrees-of-freedom of these blade tip nodes and

o i = max[u i (t)] -min[u i (t)]. (49)
• The fourth criterion evaluates the relative error on the displacement at the blade tip nodes and is defined as

e 4 = max   0, 1 - t ∆u b (t) T ∆u b (t) t u FOM b (t) T u FOM b (t)   , (50)
where

∆u b (t) = u FOM b (t) -u ROM b (t). ( 51
)
Negative values of the criteria e j (j = 1, 2, 3, 4) are set to 0 when computing the performance indicator I. Criteria e 1 and e 2 are global criteria: they allow to account for the differences in shape and in energy between the solutions obtained with the full order and the reduced order models in the entire blade. Criterion e 3 quantifies the error on the peak-to-peak amplitude of the oscillations at the blade tip. Criterion e 4 quantifies the error on the amplitude of the displacement at the blade tip nodes. Criteria e 3 and e 4 allow to account specifically for the local error at the tip. The local criteria are particularly important when considering contact simulations as contact will occur at blade tip. The weighting factors are defined such that global factors (and therefore local factors) account for 50% (w 1 = w 2 = w 3 = w 4 = 0.25).

The computation time has voluntarily not been included in the performance indicator. As shown further, the quantity of interest, i.e. the online computation time required to integrate in time the equation of motion and reproject the response in the full order space, is mainly related to the size r of the reduction basis and nearly independent from the reduction method used. The computation time is the subject of a dedicated investigation in section 4.3.2.

Relevance of the modal derivatives selection criterion

This section aims at validating the modal derivatives selection criterion introduced in section 3.2.3. A reduction basis composed of 3 constraint modes, 4 fixed interface modes and the 10 modal derivatives with the highest weights is built (see blue bars in Fig. 8). The reduced model is used to compute the response of the structure to the training excitation defined above, i.e. a harmonic excitation of amplitude A e = 300 N and pulsation ω e = 2, 000 rad/s applied at node LE in the e y direction. The transient response is shown in Fig. 9(a) and compared with the reference nonlinear solution. Qualitatively, a good match is observed between both signals. The reduction basis allows to capture the nonlinear features of the response. Quantitatively, a performance indicator I = 0.93 is obtained. The solution is also compared to the solution obtained when all the 28 modal derivatives of the 7 Craig-Bampton modes are included in the reduction basis. The addition of the remaining 18 modal derivatives in the reduction basis does not significantly enhance the performance of the reduction method. Indeed, the performance indicator is only slightly increased to I = 0.94. Lastly, the response obtained with a reduction basis composed of 10 randomly chosen modal derivatives is also represented in the figure. The performance indicator drops down to I = 0.76. The steady-state responses computed with the different reduction bases are shown in Fig. 9(b) and compared with the reference nonlinear solution. The different reduced models provide very accurate predictions and the initial errors do not build up during the simulation.

These results confirm that it is not necessary to include all modal derivatives in the reduction basis to have good performance. It also shows that the modal derivatives have to be carefully selected and underlines the relevance of the proposed selection criterion. This is particularly true when trying to accurately capture the transient response. This selection criterion is therefore a powerful tool that allows to reduce the size of the reduction basis while ensuring an accurate prediction of the structure dynamic response.

Comparison of reduction methods 4.3.1 Accuracy

The accuracy of the reduction methods is first assessed in the specific case where the external excitation is defined in the same way as the training excitation used in the offline stage for the construction of the POD reduction basis, for the POD filtering of the nonlinear forces in the Craig-Bampton method and for the selection of the modal derivatives. The harmonic excitation is characterized by an amplitude A e = 300 N, a pulsation ω e = 2, 000 rad/s and is applied at N e = 1 node (node LE). Although the training excitation is limited to a single period, the external excitation is here extended to ten periods in order to assess the capability of the reduced order models to accurately predict nonlinear transient phenomena.

The response of the structure is computed for each reduction method, for different sizes of the reduction bases, and the performance indicator is evaluated for each case. Figure 10 shows the evolution of I for each reduction method with the number r of vectors in the reduction basis. The reduction bases built by POD clearly outperform the other reduction methods in this particular case. Very good results are obtained with a small reduction basis size; the performance indicator takes values very close to 1 when r is larger than 10. The Craig-Bampton basis with POD filtering and the modal derivative approach give slightly less accurate results. Smaller reduction bases are required with the Craig-Bampton method to provide satisfactory results (performance index larger than 0.9) but it is remarkable that the modal derivative approach allows to reach higher performance indicators. For large sizes of the reduction basis, the difference in performance criteria is reflected in the global criteria: reduction bases with modal derivatives allow a better representation of the global behavior of the blade dynamics. This can be highlighted by comparing the von Mises stress fields obtained with the different approaches. Figure 11 represents the von Mises stress fields corresponding to the displacement fields predicted by the different reduction methods for a constant reduction basis size r = 25. The stress field is computed at the time t max corresponding to the second peak of the excitation. The reduction method by POD provides excellent results both locally and globally: the stress field is very close from the reference stress field computed with the full order finite element model. Regarding the two other reduction methods, this figure highlights that the solution obtained with the modal derivatives is more accurate at the scale of the entire blade than the solution obtained with the Craig-Bampton method.

Computational cost

Keeping in mind that the main objective of the reduction methods is to accelerate the simulations while ensuring a high level of accuracy, this section is dedicated to the confrontation of the three employed reduction methods in terms of computation times. It is worth reminding that the idea of reduced order modeling is to invest computational effort in the offline stage in order to reduce the computational cost of the online stage. In industrial applications, the primary objective is to minimize the computational time required by the online stage (online time) so that repeated time simulations can be performed with the same reduced model in the most efficient way. For the sake (a) of completeness, a comparative study of the computational time required by the offline stage (offline time) is also presented in this section. This allows to have an idea of the overall reduction time, which is a measure of the total effort needed to compute the response of the structure using a reduction procedure. The different computation times presented in the following correspond to single CPU sequential times on an Intel Xeon E5-2698 v4 processor (2.2 GHz). The core numerical tool is written in Python. The computation times are compared for the different reduction methods with the same reduction basis size r = 20 which leads to performance indicators larger than 0.9 for the 3 reduction methods studied.

The offline time consists in the sum of the time t 1 required to build the reduction basis, the time t 2 required to identify the nonlinear forces with the STEP and the time t 3 required to project the full order matrices onto the reduced space. The time t 3 is completely negligible with respect to t 2 . Similarly, t 1 is usually negligible with respect to t 2 , except for the POD reduction method where a full nonlinear simulation has to be performed. When POD filtering of the nonlinear forces is not considered, t 2 is proportional to r 3 , where r is the number of vectors in the reduction basis. When POD filtering is considered, the time of the full order nonlinear simulation has to be accounted for in t 2 . Table 2 compares the offline costs of the 3 studied reduction methods, for a constant basis size r = 20. As expected, the two reduction methods based on a POD require slightly more computational capabilities. It should be pointed out that the elements of tensors à and B of the same category (i.e. coefficients with equal indices, coefficients with two unequal indices and coefficients with three unequal indices) can be identified independently.

For instance, according to Eq. ( 34), the vectors Ã11 and B111 can be identified independently from the vectors Ã22 and B222 . Parallel computing can therefore be used to reduce computation times. In practice, the nonlinear force static evaluations performed on the full order finite element method can be done in parallel to build the tensors à and B, allowing to divide the computation time by the number of CPU cores or the number of software licences available if a commercial finite element software is used. The online time is mainly dominated by the time corresponding to the numerical integration of the equation of motion. Table 3 compares the online time costs of the different reduction methods (with a reduction basis size r = 20) to each other and to the time taken by the nonlinear solver of SAMCEF to numerically integrate the full order equation of motion. These computation times are measured in the framework of the simulations performed in section 4.3.1, i.e. simulations of 0.03 s with a time step of integration h = 10 -4 s. The online time is directly related to the size r of the reduction basis, and nearly independent from the reduction method used. The different reduction methods allow to reduce the computation time by a factor 2,400 with respect to the nonlinear reference model. However, as shown in the previous section, reduction bases of the same size r do not provide the same performance for the different reduction methods. For instance, reduced models built by POD allow to reach a performance indicator of 0.9 for r = 10. For such a reduction basis size, the online computation time drops down to 2 s. A compromise therefore has to be found between accuracy of the results and computation time.

Online time [s]

POD 5

CB with POD correction 5

CB and MD 5

Reference solution 12,000

Table 3. Comparison of the online computation times for the 3 reduction methods (r = 20) and the high fidelity model.

Versatility

The performance indicator I is here used to assess the robustness of each reduction method with respect to a variation of the external excitation in terms of pulsation, amplitude or spatial shape.

Pulsation and amplitude variation

The reduction bases used in Fig. 10 are reused to compute the response of the structure to a harmonic excitation of amplitude A e = 300 N and pulsation ω e = 1, 000 rad/s. The performance indicator is computed for the different methods, for each reduction basis size. The results are presented in Fig. 12.

The performances of the different reduction methods do not seem to be affected when another excitation pulsation is considered, provided the pulsation is sufficiently close to the training pulsation so that the external excitation does not significantly excite other modes of the structure. The relative performance of the methods is also conserved: the POD method provides the best results. The modal derivative method provides better results than the Craig-Bampton method with POD filtering when large sizes of the reduction basis are considered. The study of the robustness of the methods to change in the amplitude of the external excitation leads to similar conclusions.

Spatial shape variation

The robustness of the methods with respect to variation of the spatial distribution of the external excitation is then assessed. The reduction bases used in Fig. 10 and in Fig. 12 are used to compute the response of the structure to a harmonic excitation of amplitude A e = 300 N and pulsation ω e = 2, 000 rad/s distributed on the N e = 8 blade tip nodes represented in Fig. 2. The results are presented in Fig. 13. The performance of the POD method clearly drops when a different spatial distribution of the excitation is considered. The performance of the Craig-Bampton method with POD filtering is also slightly decreased. The performance of the modal derivative method remains approximately the same. A wider numerical study, considering different excitations and different numbers of physical degrees-of-freedom in the reduced space, confirms these conclusions (not shown here for the sake of conciseness). The modal derivative approach appears therefore as the most versatile reduction method among those analyzed here. Results plotted in Fig. 13 highlight the lack of versatility of the POD reduction method when the external excitation differs from the excitation used as a training to build the reduction basis. Based on this observation, a more systematic study of the influence of the shape of the training excitation has been performed. In practice, different POD bases have been built, based on different training simulations. The POD bases are composed of 3 Guyan modes and 10 POD modes obtained by SVD of the response of the structure to harmonic external excitations distributed on N t = 1, 2, . . . , 8 blade tip nodes. The eight reduction bases have then been used to compute the response of the structure to different harmonic external excitations distributed on N e = 1, 2, . . . , 8 blade tip nodes and the value of the performance indicator has been computed for each combination of POD basis and external excitation. Figure 14 represents how I evolves when the spatial distribution of the POD training excitation and the spatial distribution of the external excitation vary. As expected, the highest values of the performance indicator are obtained on the main diagonal of the figure, i.e. when the POD basis is used to compute the response of the structure to an external excitation whose spatial shape corresponds to the spatial shape of the training excitation. When the spatial shape of the external excitation differs from the training excitation applied, the performance indicator decreases. This figure allows to provide guidelines for the choice of the training excitation to get the most robust reduction method. As a general rule, the training excitation should be as close as possible to the external excitations for which the reduced model will be used. In this particular case, choosing N t = 4 seems to be a good compromise as it leads to performance indicator values higher than 0.95 when N e varies between 1 and 8. It is also interesting to notice that the method is more robust for N t = 8 than for N t = 1. In particular, using a reduction basis with N t = 8 to compute the response of the structure to an external excitation with N e = 1 provides better results than using a reduction basis with N t = 1 to compute the response of the structure to an external excitation with N e = 8.

Partial conclusion

The study of the NASA rotor 37 blade subjected to a harmonic loading and accounting for geometric nonlinearities highlights the strengths and weaknesses of the different reduction methods investigated. The model reduction by POD gives excellent results, provided that the shape of the excitation is a priori known. The Craig-Bampton method with POD filtering of the nonlinear forces provides good results locally at the blade tip, but fails to provide accurate stress fields which indicates an inaccurate representation of the displacement fields far from the boundary. The modal derivative approach provides good results, both locally and globally, and is robust with respect to changes of the external excitation. Larger sizes of reduction bases are however required to capture the nonlinear dynamics of the structure compared to the other reduction methods. In the following, only the POD and the modal derivative approaches are considered to study the blade subjected to contact interactions. The Craig-Bampton method with POD filtering of the nonlinear forces is not considered in the subsequent studies as it does not allow to capture the global dynamics of the blade.

5 Analysis with geometric nonlinearities and contact

Contact scenario

The NASA rotor 37 blade is here assumed to rotate at a constant speed Ω about the e z axis. The blade is surrounded by a casing assumed to be infinitely rigid, i.e. insensitive to blade contact events, see Fig. 15. The casing follows the shape of the blade tip, so that the clearance c 0 between the blade and the casing is uniform from the leading edge (node LE) to the trailing edge (node T E). The degrees-of-freedom of the 8 boundary nodes colored in Fig. 2 distributed between the nodes LE and T E are kept in the reduced space and define the contact interface. Possible penetrations and contact forces will be computed at these 8 nodes only. At the beginning of the simulation, the casing is perfectly circular to avoid any initial penetration, with an operating clearance between the blade tip nodes and the circular casing c 0 . In order to initiate contact, the casing is then progressively deformed in the radial direction until reaching the final distortion

f (θ) = h b exp - θ -θ/π π -π/2 0.15 2 , ( 52
)
where θ is the angular coordinate varying between 0 and 2π. This distortion corresponds to an ovalization of the casing with the creation of two diametrically opposed bumps of height h b corresponding to privileged contact areas (see Fig. 16). The gap between a given blade tip node and the casing is therefore computed through

g(θ b , t) = c 0 -(1 -e -αt )f (θ b ) -u r (t), (53)
where θ b is the current angular position of the blade tip node, u r is the radial displacement of the considered node and α is defined such that the ovalization reaches 99% of the maximal distortion at 20% of the total simulation time. In the following, the clearance c 0 is set to 0.356 mm, in agreement with the clearance of the original rotor 37 [START_REF] Moore | Performance of Single-Stage Axial-Flow Transonic Compressor With Rotor and Stator Aspect Ratios of 1.19 and 1.26[END_REF], and the height of bumps is defined as twice the clearance, h b = 0.712 mm. As detailed in section 3.4, contact is numerically handled with Lagrange multipliers and the equation of motion is integrated forward in time using an explicit central difference time integration scheme. Friction is modelled using a Coulomb law, with a friction coefficient µ = 0.15.

Preliminary studies

Contrary to the test case described in section 4 where the blade is subjected to a known harmonic excitation, it is here not possible to obtain a reference solution with the full order finite element model as the associated numerical simulations would be computationally too demanding. Performance indicators, such as the one defined by Eq. ( 43), can therefore not be directly used to assess the accuracy of the reduced model in the specific case of contact simulations. In order to give confidence in the numerical models, different preliminary studies are therefore required. Preliminary computations are carried out for a single angular speed Ω = 1, 450 rad/s for which the influence of geometric nonlinearities is observed before time and space convergence of each considered reduction method is assessed.

The responses of the structure at node LE at the angular speed Ω = 1, 450 rad/s are compared in Fig. 17 when geometric nonlinearities are taken into account and when they are neglected. The responses have been obtained with a reduced model built with the modal derivative approach, considering r b = 24 constraint modes (corresponding to the 8 boundary nodes kept in the reduced space), r c = 15 fixed interface linear modes and r MD = 40 modal derivatives. This reduction basis is significantly larger than the reduction bases used in section 4 because contact events may excite high-frequency modes. The choice of the reduction parameters is justified in section 5.2.1. During the simulation, the blade is first guided by the first bump on the casing, then oscillates freely until getting in contact with the second bump, and so on. The comparison of the steady-state responses highlights significant differences in amplitude and in frequency content between both cases and motivates the extension of existing reduction methods to account for geometric nonlinearities in contact simulations.

Time and space convergence

Convergence of the numerical procedure is first assessed with respect to the time step h of the time integration procedure. While the use of an explicit time integration scheme implies small time steps, no iterative procedure is required for contact treatment at each time step. In the present case, time steps smaller than 10 -7 s are required to ensure stable computations due to the high eigenfrequencies of the blade. Figure 18 shows the time responses obtained with the modal derivative reduction approach when considering the values h = 10 -7 s and h = 10 -8 s for the time step of integration. It confirms that the value h = 10 -7 s is sufficiently small to accurately predict the displacement of the blade. Similar conclusions hold when considering reduced order models built by POD.

The asymptotic convergence of the numerical procedure is then assessed with respect to the reduction parameters, in particular with respect to the number of modes in the reduction basis. Because contact events may yield the participation of potentially high-frequency vibration modes, a larger number of modes have to be included in the reduction basis compared to the previous study case with a harmonic excitation. Figure 19 compares the responses obtained with two reduced order models built with the modal derivative approach. The first reduction basis is composed of 24 constraint modes, 15 internal modes and 40 modal derivatives. The second reduction basis is composed of 24 constraint modes, 20 internal modes and 45 modal derivatives. The good match between the solutions, both in the transient and steady-state regimes, confirms that this reduction basis size is appropriate. When building reduced order models by POD, 20 POD modes are required to reach the asymptotic convergence of the results.

Comparison of the reduction methods

The time responses predicted with the POD and modal derivatives reduction methods are compared for the time and reduction parameters identified in the previous section. The analyses conducted in section 4 showed that the training excitation used to build the POD reduction basis has to be defined as close as possible to the actual excitation. In the case of contact studies, the excitation is not known a priori.

As a first step, the training excitation is therefore defined as a harmonic excitation in the radial direction. The response obtained with this reduced model is compared with the response obtained with the reduced model built with the modal derivative approach in Fig. 20. The responses are in good agreement during the first time steps of the simulation, but the accumulation of errors leads to very different time signals in the steady-state regime.

As a second step, the training excitation of the POD method is defined as the contact excitation predicted with the modal derivative reduced order model during the first cycle of revolution. The radial efforts at the leading edge, at mid-chord and at the trailing edge of the blade tip are represented in Fig. 21. The response obtained with this reduced model is compared with the response obtained with the reduced model built with the modal derivative approach in Fig. 22. An almost perfect agreement between the responses, both in the transient and steady-state regimes, is observed. Several conclusions can be drawn. First, the difference between the responses predicted with the POD method in Fig. 20 and in Fig. 22 confirms the extreme numerical sensitivity of the POD reduction method to the training excitation. In practice, it is necessary to know a priori the excitation to obtain good results with the POD reduction method, which is not possible in the case of contact simulations. The reduction method by POD is therefore considered not suitable for the present application. Second, this comparison gives confidence in the modal derivative approach to predict contact responses of geometrically nonlinear structures as very close results have been obtained by considering two different reduction bases built with two different approaches. The reduced model built with the modal derivative approach is therefore selected to perform more detailed studies.

Numerical investigations

The modal derivatives-based approach is used to build the reduced order model of the rotor 37 blade. Based on this model, this section details several numerical investigations, focusing on local and global quantities of interest. Based on the proposed methodology, the first analysis of the influence of geometric nonlinearities for the detection of critical interaction speeds due to contact events is presented. Radial effort computed with the modal derivative approach at the leading edge ( ), at mid-chord ( ) and at the trailing edge ( ) of the blade tip at Ω = 1, 450 rad/s with geometric nonlinearities.

Blade tip radial displacement

The first quantity of interest when analyzing contact simulations is the blade tip radial displacement as it is directly related to the occurrence of contact events. A wide angular speed range is here studied: simulations at constant angular speed Ω between 750 rad/s and 2,200 rad/s are performed (with a step of 1 rad/s). For each angular speed, 200 cycles are simulated. At each angular speed, the overall displacement level and the frequency content of the radial displacement are computed on the second half of the time signal to get rid of transient effects. The root-mean-square level of the radial displacement node LE is shown in Fig. 23 as a function of the angular speed in the cases where geometric nonlinearities are neglected and where they are considered. The frequency content of the radial displacement at node LE is then computed at each angular speed. The corresponding interaction maps are plotted in Fig. 24(a) in the case where geometric nonlinearities are not considered and in Fig. 24(b) in the case where geometric nonlinearities are taken into account. Both maps have the same color scale from blue (lower amplitudes) to red (higher amplitudes). As a first observation, it can be noted that geometric nonlinearities reduce globally the amplitude of vibrations on the whole angular speed range. This can clearly be seen in the RMS level figure (Fig. 23). In this figure, it is remarkable that both RMS signals look fundamentally different while the blade does not a priori seem to undergo large displacement (radial displacement of the order of 1 mm). These differences in RMS levels show that geometric nonlinearities may have a strong influence on blades contact dynamics, even for blades that are not particularly slender, such as compressor blades.

Both interaction maps show lines of higher amplitudes on even engine order lines, which is consistent with the fact that the casing ovalization leads to an excitation occurring twice per revolution. In the case where geometric nonlinearities are taken into account, richer frequency content are observed, for instance with contributions of odd harmonics visible around 1,000 rad/s. Non-integer harmonics of the angular speed are also visible around 2,100 rad/s in Fig. 24 The interaction maps and the RMS level figures show the existence of strong contact interactions between the first vibration modes of the blade and even engine orders. These interactions are characterized by high amplitudes of vibration followed by amplitude jumps characteristic of nonlinear systems where the amplitude of the displacement suddenly drops for a small change of the angular speed. For the case without geometric nonlinearities, the interaction map highlights the existence of interactions between the first bending mode and the engine orders 8, 6 and 4 (denoted by 1B/e o8 , 1B/e o6 and 1B/e o4 ) at the angular speeds Ω = 872 rad/s, Ω = 1, 209 rad/s and Ω = 2, 074 rad/s. Because of contact stiffening, the angular speeds of the interactions are slightly higher than the angular speeds corresponding to the exact intersections between the first bending mode and the engine order lines, respectively Ω = 712 rad/s, Ω = 949 rad/s and Ω = 1, 424 rad/s. Geometric nonlinearities are responsible of an additional stiffening. As shown in the interaction map in Fig. 24(b), the interactions 1B/e o8 , 1B/e o6 and 1B/e o4 occur at higher angular speeds, respectively Ω = 1, 145 rad/s, Ω = 1, 682 rad/s and outside the studied angular speed range. To the best of the author's knowledge the stiffening due to the combined effect of geometric and contact nonlinearities has never be characterized on industrial mechanical structures. Several studies have however shown that the stiffening caused by large displacements can reach 10% to 40% in term of frequencies [START_REF] Georgiades | Modal Analysis of a Nonlinear Periodic Structure with Cyclic Symmetry[END_REF][START_REF] Grolet | Vibration Analysis of a Nonlinear System With Cyclic Symmetry[END_REF][START_REF] Martin | Dynamic Analysis and Reduction of a Cyclic Symmetric System Subjected to Geometric Nonlinearities[END_REF]. Some interactions with other vibration modes also lead to high amplitude displacements, such as interaction 1T/e o12 at Ω = 1, 492 rad/s and interaction 2B/e o10 at Ω = 2, 181 rad/s for the case without geometric nonlinearities and interaction 2B/e o14 at Ω = 1, 789 rad/s in the case with geometric nonlinearities.

These interaction maps also highlight the existence of more complex interactions involving several vibration modes. At Ω = 1, 859 rad/s, for instance, the geometrically linear solution can be written as stable combination of the first bending mode and the first torsion mode. In the case with geometrical nonlinearities, two stable solutions coexist between Ω = 1, 550 rad/s and Ω = 1, 682 rad/s, one of high amplitude dominating by the first bending mode and one of smaller amplitude dominating by first torsion mode.

The direct analysis of the displacement interaction maps only allows to have a global view of the interactions. The interactions could be further characterized by using frequency methods or looking at the modal components but this is not the subject of this article.

Stress analysis

From an industrial point of view, the stress field inside the blade is an important global quantity of interest as it is directly related to the structure integrity. If the maximal stress level exceeds the material yield stress, irreversible deformations may occur, leading eventually to the failure of the structure. The stress fields are here studied at the angular speed corresponding to the interaction 1B/e o8 which is visible in both interaction maps, i.e. at Ω = 872 rad/s in the case where geometric nonlinearities are neglected and at Ω = 1, 145 rad/s in the case where geometric nonlinearities are considered. These angular speeds are indicated by the symbol ♠ in Fig. 24. The von Mises stress fields computed at time t max , different for both cases, where the node LE reaches its maximal radial displacement are represented in Fig. 25. At time t max , the blade is not subjected to any contact force and is free to oscillate between two bumps.

At these interactions, the response of the structure is dominated by the first bending mode of the blade, which shows slightly different shapes when geometric nonlinearities are considered or not. The zones of maximal von Mises stresses are not located at the same place in both cases. When geometric nonlinearities are neglected, the maximal stresses are found near the root of the blade. When considering the geometric nonlinearities, the zone of maximal stresses is moved toward the center of the blade. These results suggest that accounting for geometric nonlinearities in bladed disk models could have a significant influence in term of blades design.

Quantitatively, the ratio of the maximal von Mises stresses in both cases is in line with the corresponding ratio of maximal displacements at the blade tip leading edge (about 2). This confirms that, in the studied case, geometric nonlinearities weakens the interactions, not only locally at the blade tip, but also globally in the whole blade. The analysis of the stress field at the other interaction angular speeds confirms that the severity of the interaction in terms of displacement amplitudes is reflected in the maximal values of the von Mises stress. Neglecting the geometric nonlinearities leads to an overestimation of the stress level in the blade.

Conclusion

This paper presents the set up of a new methodology to predict the dynamics of geometrically nonlinear structures subjected to contact interactions in a numerically efficient way through the use of reduced order models. In order to develop this methodology, the use of different reduction bases has been investigated following respectively the proper orthogonal decomposition approach, the Craig-Bampton approach and the modal derivative approach. The different reduction methods have been adapted if required and used to study the dynamics of a compressor blade from the NASA rotor 37 in the case of a harmonic excitation, representative of an aerodynamic loading, and in the case of contact excitation. These studies allow to highlight the strengths and weaknesses of the different methods. The reduction method by POD proves to be very sensitive to the spatial shape of the excitation. In contact simulations, where the exact excitation is not known a priori, POD does not provide accurate results. The Craig-Bampton reduction method does not allow to capture globally the dynamics of the blade, even in the simple case of a harmonic excitation. The modal derivative approach allows to build robust reduced order models. In order to limit the size of the reduced system, a modal derivative selection criterion has been proposed.

The final methodology proposed here to predict the nonlinear and nonsmooth dynamics of the blade in a computationally efficient way is based on the projection of the high fidelity model onto a reduction basis composed of Craig-Bampton modes and a selection of their modal derivatives. The nonlinear internal forces due to large displacements are evaluated with the stiffness evaluation procedure. Contact is numerically handled with Lagrange multipliers and the equation of motion is integrated forward in time using an explicit central difference time integration scheme. This methodology has been applied to study in more details the contact interactions of the rotor 37 blade with a surrounding rigid casing. On the one hand, the numerical studies have shown that the reduction procedure developed allows to extract quantities of interest that are relevant for both researchers and industrial designers. In particular, it can be used to efficiently identify the critical angular speeds at which high-amplitude interactions occur and to asses the consequences of the interactions on the structure integrity, for instance in terms of stress fields inside the blade. On the other hand, the comparison of the results obtained without and with geometric nonlinearities shows that neglecting the geometric nonlinearities in such studies can lead to wrong estimation of the critical angular speeds of the structure and to bad evaluation of the severity of the interactions in terms of displacement amplitudes and von Mises stress fields.

Figure 1 . 1 Figure 2 .

 112 Figure 1. Picture of NASA rotor 37 [43].

Figure 3 .

 3 Figure 3. First three free vibration modes of the NASA rotor 37 blade: (a) mode 1 (1B), (b) mode 2 (1T) and (c) mode 3 (2B).

  0

Figure 5 .

 5 Figure 5. Linear ( ) and nonlinear ( ) responses of the blade at node LE in the ex direction (frequency-domain).

Figure 6 .

 6 Figure 6. Offline and online stages in model order reduction.

Figure 7 .

 7 Figure 7. Comparison of the nonlinear responses of the blade at node LE in the ex direction obtained with the high fidelity model ( ) and with the Craig-Bampton reduced model without POD correction ( ).

Figure 8 .

 8 Figure 8. Weights Wij associated to each modal derivative (10 highest weights ( )).

Figure 9 .

 9 Figure 9. Comparison of the reference nonlinear response of the structure ( ) with the responses obtained with all modal derivatives ( ), with the best 10 modal derivatives ( ) and with 10 randomly chosen modal derivatives ( ): (a) transient response and (b) steady-state response.

Figure 10 .

 10 Figure 10. Evolution of the performance indicator for the different reduction methods (POD (•), Craig-Bampton with POD filtering ( ), modal derivatives ( )) as a function of the size of the reduction basis (Ae = 300 N, ωe = 2, 000 rad/s, Ne = 1).

1 Figure 11 .

 111 Figure 11. Comparison of the von Mises stress fields at tmax obtained with the full order finite element model, the reduced model obtained with Craig-Bampton and POD filtering and the reduced model obtained with the modal derivative approach (r = 25): (a) reference solution, (b) POD, (c) CB with POD filtering and (d) modal derivatives.

Figure 12 .

 12 Figure 12. Evolution of the performance indicator for the different reduction methods (POD (•), Craig-Bampton with POD filtering ( ), modal derivatives ( )) as a function of the size of the reduction basis (Ae = 300 N, ωe = 1, 000 rad/s, Ne = 1).

Figure 13 .

 13 Figure 13. Evolution of the performance indicator for the different reduction methods (POD (•), Craig-Bampton with POD filtering ( ), modal derivatives ( )) as a function of the size of the reduction basis (Ae = 300 N, ωe = 2, 000 rad/s, Ne = 7).

Figure 14 .

 14 Figure 14. Performance indicator as a function of the number of nodes excited in the POD training (Nt) and excited by the external excitation (Ne). Results obtained with POD bases composed of 3 Guyan modes and 10 POD modes.

1 Figure 15 .

 115 Figure 15. Cut view of the NASA rotor 37 blade surrounded by the casing and clamped on the disk.

1 Figure 16 .

 116 Figure 16. Progressive ovalization of the casing (not to scale).

Figure 17 .Figure 18 .

 1718 Figure 17. Comparison of the responses of the blade at node LE in the radial direction to a contact excitation at Ω = 1, 450 rad/s with geometric nonlinearities ( ) and without geometric nonlinearities ( ): (a) transient response and (b) steady-state response.

Figure 19 .

 19 Figure 19. Influence of the reduction parameters on the response of the blade at node LE in the radial direction to a contact excitation at Ω = 1, 450 rad/s with geometric nonlinearities (reduction bases composed of 24 contraint modes, 15 internal modes and 40 modal derivatives ( ); 24 contraint modes, 20 internal modes and 45 modal derivatives ( )): (a) transient response and (b) steady-state response.

Figure 20 .Figure 21 .

 2021 Figure 20. Comparison of the response of the blade at node LE in the radial direction to a contact excitation at Ω = 1, 450 rad/s with geometric nonlinearities (reduction bases built with the modal derivative approach ( ) and by POD approach with a harmonic training excitation ( )): (a) transient response and (b) steady-state response.

Figure 22 .2Figure 23 .

 2223 Figure 22. Comparison of the response of the blade at node LE in the radial direction to a contact excitation at Ω = 1, 450 rad/s with geometric nonlinearities (reduction bases built with the modal derivative approach () and by POD with the contact excitation computed with the modal derivative approach ()): (a) transient response and (b) steady-state response.

  (a) and around 1,550 rad/s in Fig. 24(b), which can be related to the nonlinear nature of the dynamics. Delhez et al.

Figure 24 .

 24 Figure 24. Interaction maps of the radial displacement at node LE (with eigenfrequencies ( ) and engine order lines ( )): (a) without geometric nonlinearities and (b) with geometric nonlinearities.

Figure 25 .

 25 Figure 25. Comparison of the von Mises stress fields at tmax obtained with the full order finite element model, the reduced model obtained with Craig-Bampton and POD filtering and the reduced model obtained with the modal derivative approach (r = 25): (a) without geometric nonlinearities and (b) with geometric nonlinearities.

Table 1 .

 1 Characteristics of the reduced bases.

Table 2 .

 2 Comparison of the offline computation times for the 3 reduction methods (r = 20).

Delhez et al.

Acknowledgments

E. Delhez is supported by the Fonds de la Recherche Scientifique (F.R.S.-FNRS, Belgium) which is gratefully acknowledged. This research was also undertaken thanks to funding from the Canada Research Chairs Program. Delhez et al.