Nicolas Di Palma
email: nicolas.di-palma@ec-lyon.fr

Benjamin Chouvion

Fabrice Thouverez

Parametric Study on Internal Resonances for a Simplified Nonlinear Blade Model

Keywords:

This article proposes to determine and identify the physical contributions that drive the internal resonance between the bending and torsion modes of a blade. The blade is modeled as a nonlinear beam including warping, twist and pre-bending, and a discretization scheme based on a Rayleigh-Ritz approximation is performed. The modal properties of the model are compared and validated using beam theory and finite element modeling. Several configurations of the model allowing the occurrence of an internal resonance between the first torsional and the second bending modes are considered. The frequency responses near the second bending mode are obtained through the Harmonic Balance Method coupled with a branch switching algorithm to locate and track bifurcated branches related to internal resonances. Attention is paid to the identification of the geometrical nonlinear terms that enable the appearance of an internal resonance. With the help of a multiple scale analysis, it is shown that it is essentially driven by a small set of coupling parameters involving bending and torsion displacements of the commensurate modes. It is also observed in a rotating situation thanks to the hardening of the bending mode.

Introduction

In response to the current environmental and economic challenges, aircraft engine manufacturers are constantly looking toward new efficient designs. One field of development concerns the conception of engines with a larger diameter. A direct consequence of this is an increase in the size of engine components. This is true of fan blades, for which slender geometries are derived. Such a design facilitates the emergence of geometrical nonlinearities inducing large deformations and displacements in the structure. This nonlinear behavior can manifest itself in different aspects.

One of them is referred to in the literature as internal resonance. It expresses itself through modal coupling whose respective frequencies are commensurate [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Due to the amplitude-frequency dependency in nonlinear systems, internal resonance can stem from systems that do not necessarily have linear commensurate frequencies [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]. In practical terms, it manifests itself by a multi-modal response and leads to unexpected large vibration levels. Hence, it is undesirable and represents a serious issue in the conception phase as it can reduce components' lifetime and increase fatigue risks.

Internal resonances have been broadly studied in the literature. Due to its simple implementation, the analytical method of multiple scales has been intensively employed to study internally resonant systems having a small number of degrees of freedom (dofs). In their reference book [START_REF] Nayfeh | Nonlinear oscillations[END_REF], Nayfeh and Mook addressed a large number of internally resonant cases in academic models with cubic and quadratic nonlinearities. Nayfeh et al. [START_REF] Nayfeh | Nonlinear coupling of pitch and roll modes in ship motions[END_REF] observed a saturation phenomenon for a 2-dofs system when studying a 1:2 internal resonance with multiple scales analysis. Additionally, cascades of internal resonances involving three modes have also been investigated [START_REF] Monteil | Nonlinear forced vibrations of thin structures with tuned eigenfrequencies: the cases of 1: 2: 4 and 1: 2: 2 internal resonances[END_REF] using this approach. Though it furnishes interesting insights, multiple scales analysis remains confined to weak nonlinearities and restricted to a small number of dofs. Methods relying on nonlinear normal modes (NNMs) have also been proposed to study internal resonances. Works have suggested extending the notion of invariant manifold to a multi-mode invariant manifold by wisely selecting system state variables involved in the internal resonance [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF][START_REF] Pesheck | Nonlinear modal analysis of structural systems using multi-mode invariant manifolds[END_REF]. Other approaches including the shooting method have been used to compute the NNMs of systems involving internal resonances [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF]. It has been shown that internal resonances are manifested by tongues on frequency energy plots when they reach a sufficient level of energy [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]. They have also been computed and studied using Galerkin methods such as the Harmonic Balance Method (HBM). Lewandowski [START_REF] Lewandowski | Non-linear free vibrations of beams by the finite element and continuation methods[END_REF][START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structures-part 2: numerical strategy and examples[END_REF] used the HBM to compute free and forced responses of beams in presence of internal resonances. In [START_REF] Lau | Nonlinear vibration of thin elastic plates, part 2: internal resonance by amplitude-incremental finite element[END_REF], internal resonances were observed in the study of large-amplitude periodic vibrations of thin elastic plates. Additionally, Ribeiro and Petyt proposed to combine the finite element approach with the HBM and observed internal resonances for a uniform slender beam [START_REF] Ribeiro | Non-linear vibration of beams with internal resonance by the hierarchical finite-element method[END_REF] and isotropic plates [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF]. In [START_REF] Stoykov | Nonlinear free vibrations of beams in space due to internal resonance[END_REF], the HBM was employed on beam equations based on Timoshenko's and Saint Venant's formulations and it was demonstrated that internal resonances lead to bending-torsion coupling.

Although internal resonances have been extensively studied, efforts still need to be made to fully understand the phenomenon, especially in the field of turbomachinery (see e.g. [START_REF] Palma | Nonlinear harmonic analysis of a blade model subjected to large geometrical deflection and internal resonance[END_REF]). More specifically, to the authors' knowledge, no in-depth studies focusing on the identification of the physical quantities responsible for the appearance of internal resonance have been carried out to date. These observations are the reason behind the present work, which proposes to perform a parametric study on internal resonances for a simplified blade model.

To address this issue, an accurate model of the blade restituting the nonlinear effects needs to be derived. A finite element procedure would give the most accurate blade model. However, it requires dealing with large dimension variables. However, a phenomenological model would represent a good compromise between accuracy and number of dofs. In addition, it enables identification of the effects of each parameter independently. Many efforts have been made to obtain a nonlinear blade model that is simplified but with a sufficient level of accuracy. Houbolt and Brooks [START_REF] Houbolt | Differential equations of motion for combined flapwise bending, chordwise bending, and torsion of twisted nonuniform rotor blades[END_REF] derived the equations of motion of nonuniform twisted blades including bending-torsion coupling. However, they neglected the nonlinear terms and only retained linear contributions. Hodges and Dowell [START_REF] Hodges | Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[END_REF] improved on that work by retaining the dominant nonlinear terms using an ordering scheme. Further contributions have been proposed to improve the nonlinear equations of motion for blades [START_REF] Kaza | Nonlinear curvature expressions for combined flapwise bending, chordwise bending, torsion and extension of twisted rotor blades[END_REF][START_REF] Rosen | Nonlinear equations of equilibrium for elastic helicopter or wind turbine blades undergoing moderate deformation[END_REF][START_REF] Crespodasilva | Flag-lag-torsional dynamics or extensional and inextensional rotor blades in hover and in forward flight[END_REF][START_REF] Subrahmanyam | Nonlinear vibration and stability of rotating, pretwisted, preconed blades including coriolis effects[END_REF]. Recently, Riziotis et al. [START_REF] Riziotis | Aeroelastic analysis of pre-curved rotor blades[END_REF] considered a pre-curved twisted rotor blade in which the pre-bending was modeled by a change in the elastic axis of the blade. Their analysis relied on two approaches: a multi-body approach and a method similar to the one employed in [START_REF] Hodges | Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[END_REF]. However, warping of the blade cross-section was voluntarily omitted in this work to simplify the derivation of the equations. In this paper, extending the work presented in [START_REF] Hodges | Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[END_REF], we derive the equations of motion of a simplified nonlinear blade model including pre-twist, pre-bending, and warping. Pre-bending enables the slender shape of a fan blade to be reproduced. The equations of the model are expanded up to the second order and the frequency response is computed thanks to the HBM. This choice is suitable as it allows us to obtain multi-mode responses in case of commensurability.

Besides, the HBM can be coupled to a branch switching algorithm that finds and tracks bifurcated branches related to internal resonances.

The originality of our work is that we perform a parametric study on different configurations of the blade model to derive the physical quantities driving the internal resonance between bending and torsion modes. Special attention will be paid to the influence of cross-section shape, pre-bending, and rotational speed on the appearance of internal resonance.

The paper is organized as follows. Section 2 introduces the twisted and pre-bent simplified blade model and establishes its differential equations of motion employing a Rayleigh-Ritz discretization procedure. A validation scheme is included to assess the precision of the modeling. Then, section 3 recalls the harmonic balance formulation and the branch switching algorithm used to track branches of solutions emerging from internal resonances. Finally, section 4 details the results of the parametric analysis associated with the internal resonance and based on the sensitivity of the geometrical contributions.

Construction of a second order pre-bent twisted blade model

The blade model presented in this section is derived from the work of Hodges and Dowell [START_REF] Hodges | Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[END_REF], to which an initial pre-bending is added following the approach of Riziotis et al. [START_REF] Riziotis | Aeroelastic analysis of pre-curved rotor blades[END_REF]. The model consists of a nonlinear pre-bent and twisted rotor blade including warping of the cross-sections and is expanded up to the second order. The kinematics of the model is first described. The kinetic and potential energies are then calculated. Finally, the equations of motion are derived using a Rayleigh-Ritz discretization and Lagrange's equations.

Local and global coordinate systems

Let (O,x,y,z) be a global coordinate system rotating around the z-axis at constant angular speed Ω. A blade with an arbitrary cross-section and clamped at origin O is considered. It is rotating at the same angular speed Ω around the z-axis. The centerline of the blade is described by its elastic axis ( ) which defines an initial pre-bending in the (O,x,z) plane (see figure 1). The initial twist of the blade is defined by the function ( ) (see figure 2). Blade deformations are characterized by the displacements ( , ), ( , ) and ( , ) parallel to the x, y and z axes, as shown in Figure 1. These displacements are applied on the elastic axis of the undeformed blade so that a point located at ( , 0, ( )) moves to position ( + , , ( ) + ) after deformation. The deformed state is defined by a local coordinate system (O',x',y',z'). The x' axis is tangent to the deformed elastic axis at point O', while the y' and z' axes correspond to the principal axes of the cross-section. Their intersection O' is the shear center of the deformed cross-section. Coordinates on the principal axes are denoted respectively and (see Figure 2). In the following, for the sake of conciseness, the cross-section is assumed to be at least symmetric about the y' axis.

x O y z u(x, t) v(x, t) w(x, t) x' y' z' O' Ω z e (x)
The deformed state is also described by the torsion angle ( , ). Therefore, the twist angle of the blade changes from ( ) in the undeformed state to ( ) + ( , ) for the deformed blade. In the remainder of this paper, in order to be consistent with an engine configuration, we assume ( = 0) = ∕2. The complete motion of the blade is finally described by the translations ( , ), ( , ), and ( , ) of its elastic axis and the angle ( , ) about this axis. The dependency in and will be omitted in the following analytical developments. By using the previous definitions, the vector position 1 of an arbitrary point of the deformed blade can be expressed by: 1. The position of the deformed elastic axis in global coordinates: + , , + ;

y z v(x,t) w(x,t) z e (x) y',η O' z',ζ θ(x) θ(x) + φ(x,t)
2. The position relative to the deformed elastic axis in local coordinates: -( + ) + , , ;

where (.) + designates the derivative with respect to the curvilinear coordinate along the deformed blade. The term ( + ) + mentioned above refers to an axial position involving the warping function ( , ). This function is antisymmetric and only depends on the coordinates and of the cross-section. Introduced by De Saint Venant in his manuscript [START_REF] De Saint-Venant | Mémoire sur la torsion des prismes: avec des considérations sur leur flexion ainsi que sur l'équilibre intérieur des solides élastiques en général: et des formules pratiques pour le calcul de leur résistance à divers efforts s'exerçant simultanément[END_REF], the warping function authorizes the plane cross-section not to remain plane during deformation. In practical terms, the warping function is the solution of a Laplace equation coupled with a Neumann limit condition (see Appendix A).

The vector position 1 is therefore given by

1 = ⎡ ⎢ ⎢ ⎣ + + ⎤ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎣ -( + ) + ⎤ ⎥ ⎥ ⎦ (1)
where is the transformation matrix relating global (O,x,y,z) to local (O,x',y',z') coordinates. Its expression is given in Appendix B.

Ordering scheme

An ordering scheme is adopted to neglect higher-order terms in the equations. It relies on magnitude assumptions and is performed by introducing a small dimensionless parameter such that each dimensionless variable of the model can be expressed in terms of . The assumptions made by Hodges and Dowell in [START_REF] Hodges | Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[END_REF] are retained and it is assumed that the dimensionless pre-bending is of the same order as the dimensionless bending displacements. The different orders used in the following developments are recalled in Table 1.

Table 1

Orders of the dimensionless parameters. L is the blade length,

= and = ∕ O( 2 ) ∕ O( ) ∕ O( ) ∕ O( ) ∕ O( ) ∕ 2 O( 2 ) O( ) ∕ O( ) ∕ O(1) ∕ O( ) ∕ O( ) O (1)
In the following, the analytical developments are expanded up to the second order in . Some higher-order terms will be kept for numerical conditioning purposes and a better physical evaluation. They will be described in the following sections. For moderate pre-bending, note also that d.∕ d ≃ d.∕ d to the second order (see Appendix B).

Expression of the strain tensor

A general expression of the Green's strain tensor is determined by means of:

d 1 ⋅ d 1 -d 0 . d 0 = 2 d d d ⎡ ⎢ ⎢ ⎣ d d d ⎤ ⎥ ⎥ ⎦ (2)
where d 1 is the differential of the position vector 1 of an arbitrary point ( 0 is the position vector of the same point in the undeformed state). The variables d , d , and d represent increments along the deformed elastic axis and the two principal axes of the cross-section.

Assuming zero transversal stresses and and shearing stress , the substitution of (1) in (2) yields the following second-order approximated strain-displacement expressions:

= ′ + 1 2 ′ 2 + 1 2 ′ 2 + ′ ′ -′′ + 2 + 2 ′ ′ + 1 2 ′ 2
-′′ [ cos( + ) -sin( + )] -′′ + ′′ [ sin( + ) + cos( + )]

+ ′′ [ sin + cos ] (3a) = -( + ) ′ (3b) = ( -) ′ (3c)
In addition to the classical tension and bending terms, note the presence of further terms involving twists , warping 135 , and pre-bending in the axial strain expression (3a). The third-order term 2 + 2 ′2 in equation (3a) allows us to take into account tension-torsion coupling in the motion equations (see [START_REF] Houbolt | Differential equations of motion for combined flapwise bending, chordwise bending, and torsion of twisted nonuniform rotor blades[END_REF]).

Derivation of potential and kinetic energies 2.4.1. Potential energy

Considering the standard Hooke's law and assuming an isotropic and homogeneous model, the elastic strain energy of the blade is given by

= 1 2 ∫ Ω 2 + 2 + 2 dΩ (4)
where is the Young's modulus and = 2(1+ ) is the shear modulus. Substituting the expressions (3a), (3b), and (3c) into the elastic strain energy (4) yields to the second order:

= 1 2 ∫ 0 ′2 + 1 4 ′4 + 1 4 ′4 + ′ ′2 + ′ ′2 + 1 2 ′ ′ 2 + ′ ′ 2 + 2 ′ + ′2 + ′2 ′ ′ - ′′ (cos -sin ) + ′′ (sin + cos ) + ′′ cos 2 ′ + ′2 + ′2 + 2 ′ ′ + 2 ′ ′ + 1 2 ′2 2 ′ + ′2 + ′2 + 2 ′ ′ + ′ 1 2 ′′2 (
′ ′ 2 + 1 ′′2 + ′2 d (5)
where the following constants are introduced:

≡ ∬ d d ′ ≡ ∬ 2 d d ≡ ∬ d d ≡ ∬ ( -) 2 + ( + ) 2 d d * 2 ≡ ∬ ( 2 + 2 ) d d 2 ≡ ∬ ( 2 + 2 ) d d ′ ≡ ∬ 2 d d * 1 ≡ ∬ ( 2 + 2 ) 2 d d * 1 ≡ ∬ d d 1 ≡ ∬ 2 d d (6)
One may recognize some classical constants such as the second moments of area ′ and ′ and the polar moment of 140 area 2 . Some constants that depend on the warping function are also introduced. This is the case for the torsion constant , the warping rigidity 1 , and the two constants * 1 and * 2 . The constant characterizes the distance between the geometric center of the cross-section and the shear center (i.e the elastic axis).

Higher-order terms (underlined in ( 5)) must be kept in our approximation in order to better evaluate the torsion modes and their evolution due to the centrifugal effects.

The derived elastic strain energy is valid for a slender, moderately pre-bent, and twisted blade. It allows spanwise variation of the geometric properties (area, second moments of area, etc.).

Kinetic Energy

The kinetic energy of the rotating blade is defined by

= 1 2 ∫ Ω ⋅ dΩ (7)
with the density and the velocity of an arbitrary point on the deformed blade, defined by

= d d + ∧ 1 , = [0 0 Ω] (8)
Substituting equation ( 8) in ( 7) yields to the second order

= 1 2 ∫ 0 ̇ 2 + ̇ 2 + ̇ 2 + 2 1 + 2 2 ̇ 2 -2 ̇ ̇ sin + 2 ̇ ̇ cos + Ω -2 ̇ ( + cos ) + 2 ̇ ′ cos + 2 2 cos 2 + 2 1 sin 2 + 2 ̇ ′ sin 2 2 -2 1 cos sin + 2 ̇ -2 ̇ ′ cos -2 ̇ ′ sin -2 ̇ (sin + cos ) -2 ̇ ′ sin + Ω 2 2 + 2 (cos -sin ) + 2 -2 ′ (sin + cos ) + 2 2 2 1 + 1 -2 2 cos (2 ) -2 sin (2 ) -2 ′ (cos -sin ) + 2 1 2 1 -1 -2 2 cos (2 ) + 2 sin (2 ) -2 ′ (sin + cos ) d (9)
where the following constants are introduced:

≡ ∬ d d ≡ ∬ d d 2 1 ≡ ∬ 2 d d 2 2 ≡ ∬ 2 d d (10)
The linear mass of the blade is . The constants represent the moments of inertia with respect to the principal axes of the cross-section. Their sum is the polar moment of inertia. Finally, corresponds to the offset 150 between the center of mass and the shear center of the cross-section.

Again, four doubly underlined terms corresponding to higher-order terms are conserved in equation [START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structures-part 2: numerical strategy and examples[END_REF]. The first two prevent the matrix from not being definite and positive while the last two allow a better evaluation of the spin-softening matrix .

Work of external forces 155

In the following, no body forces are considered. Thus, the work of external forces s is given by the following expression:

= ∫ Γ ⋅ s dΓ (11)
with = [ , , , ]. In the remainder of this paper, only localized harmonic forces applied at the blade tip are considered. For a harmonic force imposed at = , the work of the external forces has the following form:

= ( ) cos( ) (12)
where and 0 are respectively the frequency and the amplitude of the external excitation.

Discretization through a Rayleigh Ritz approach

In this paper, a discretized form of the equations of motion is proposed. To this end, the Rayleigh-Ritz approach is employed and the previously introduced displacements are approximated according to the following expressions:

( , ) = u ∑ =1 ( ) ( ) = u u ( , ) = v ∑ =1 ( ) ( ) = v v ( , ) = w ∑ =1 ( ) ( ) = w w ( , ) = ∑ =1 ( ) ( ) = (13)
where the (respectively , and ) are kinematically admissible functions weighted by the generalized coordinates (respectively , and ). The constants u , v , w and correspond to the number of Ritz functions chosen to approximate the displacements , , and angle . In order to satisfy the clamped condition at the blade root = 0, the following polynomial Ritz functions are introduced:

( ) = ∈ 1, u ( ) = +1 ∈ 1, v ( ) = +1 ∈ 1, w ( ) = +1 ∈ 1, (14)
Noting = , the discretization procedure is provided by means of the Lagrange equations defined by:

d d  ̇ -  = 0 = 1, … , u + v + w + (15)
with  = -+ . Application of the Lagrange equations ( 15) leads to the following nonlinear differential equation of motion:

̈ + ( + ) ̇ + ( + ) + nl ( ) = Ω + ext (16)
with

nl ( ) = q ⊗ + c ⊗ ⊗ ( 17
)
and where

= 2 | | | | = , q = 1 2 3 | | | | = , c = 1 6 4 | | | | = ext = (18) = ̈ d d ̇ - , = ̇ d d ̇ - = d d ̇ - , Ω = d d ̇ - | | | | = , ̇ = , ̈ = (19)
Equation ( 16) represents the discrete nonlinear equation of motion of the pre-twisted and pre-bent blade model. It embeds Ω-dependent matrices and related respectively to the spin-softening and gyroscopic effects. The vector Ω corresponds to the centrifugal efforts. Moreover, the geometrical nonlinear effects are taken into account through matrices q and c giving respectively quadratic and cubic stiffnesses. The damping matrix introduced in equation ( 16) will be detailed later.

Validation of the model

The objective of this section is to validate the proposed model. The validation procedure is based on modal analyses of a simple straight beam and a twisted pre-bent rotating beam.

First of all, the structure considered is a non-rotating straight (no twist, no pre-bending) beam with a rectangular cross-section of area w×h. The parameters used are given in Table 2.

Table 2 Parameters of the straight beam model for the first validation step

Geometry Material

Length

L [m] Width w [m] Height h [m] [kg m -3 ] E [GPa] 1.3 0.4 0.04 1550 100
The first seven natural frequencies found with the proposed model are compared in Table 3 with the ones obtained from the Euler-Bernoulli theory and whose analytical expression can be found, for instance, in [START_REF] Blevins | Formulas for natural frequency and mode shape[END_REF]. For the proposed model, the Rayleigh-Ritz discretization is performed by considering 6 approximation functions per displacement ( u = v = w = = 6). The frequencies of both models match almost perfectly, as the relative error is much lower than one percent for the first few modes of interest.

A second validation step is carried out including twist, pre-bending, and rotation of the beam. The cross-section remains rectangular and the parameters used are indicated in Table 4.

Our proposed modeling is compared with a finite element simulation carried out on Ansys and based on the same geometrical and material properties. To ensure modal convergence, the FE model is meshed with 3D block elements considering 40 elements spanwise and 10 by 2 elements on the cross-section, leading to a total of 14,405 dofs. Each simulation will be performed by taking into account the geometric nonlinear effects. For this second validation step, the modes are computed near a prestressed state induced by the centrifugal effects. The prestressed state Ω is obtained by prescribing ( ) = Ω + ( ) in equation ( 16) and then solving the following nonlinear static equation:

( + ) Ω + nl ( Ω ) = Ω (20)
The modes are finally computed by solving the generalized eigenvalue problem

Λ 2 + Λ + t = (21)
with t the tangent stiffness matrix defined for the prestressed state by:

t = + + nl Ω (22)
The evolution of the natural frequencies, found with both methodologies, with the rotation is depicted in Figure 3.

A satisfactory agreement can be observed between the two models, especially for the first modes. Additionally, the trends of the different frequencies' evolutions are also respected. This is particularly the case for the second (second bending) and third (first torsional) modes, for which one can observe a veering zone near 4000 RPM.

To complete these observations, Figures 4 and5 illustrate the mode shapes at 0 RPM and 5000 RPM of the second bending and first torsional modes for both the proposed model and the Ansys model. At both rotating speeds, a very good agreement can be observed for each of the two modes. The deformed shapes obtained at 5000 RPM are even more encouraging as they show that the proposed model can accurately catch the coupling between the modes occurring in this veering area.

These validation steps instill confidence in the ability of the proposed model to represent both the linear and the nonlinear dynamic behavior of blade-type structures.

Numerical Resolution

This section recalls the main steps of the Harmonic Balance Method (HBM) and describes the branch-switching algorithm used to track and follow the bifurcated branches of solutions. A multiple scale analysis is also employed on a small phenomenological model to study the nonlinear coefficients involved in the internal resonance phenomenon. Numerous techniques abound in the literature for studying nonlinear systems like that exhibited by equation ( 16). When searching for steady-state periodic solutions, the HBM is a valuable tool [START_REF] Krack | Harmonic balance for nonlinear vibration problems[END_REF]. It has proven its great efficiency in dealing with nonlinear systems including friction nonlinearities [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF][START_REF] Chouvion | Vibration analysis of beam structures with localized nonlinearities by a wave approach[END_REF][START_REF] Firrone | The effect of underplatform dampers on the forced response of bladed disks by a coupled static/dynamic harmonic balance method[END_REF] or geometric nonlinearities [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF][START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] as well as non-smooth mechanical systems [START_REF] Kojtych | Experimental and numerical modelling of a one-degree of freedom non-smooth mechanical system[END_REF]. It has also shown great potential to study internal resonance [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF][START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF]. The HBM is a Galerkin method whose cornerstone is to assume the solution of ( 16) as a truncated Fourier series defined by:

( ) = 0 + h ∑ =1 cos( ) + sin( ) (23)
with corresponding to the excitation frequency of the system. Substituting ( 23) in ( 16), the Galerkin procedure [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF] consists in projecting the resulting expression on the Fourier basis using the following scalar product:

⟨ ( ), ( )⟩ = ∫ 2 0 ( ) ( )d (24)
This last step yields the following nonlinear algebraic system:

( ̂ , ) ∶= ( ) ̂ + ̂ nl ( ̂ ) -̂ ext = ( 25
)
where ̂ is the new vector of unknowns filled by the harmonic coefficients 0 , 1 , 1 … , and is the dynamic stiffness matrix defined by

= ⎡ ⎢ ⎢ ⎢ ⎣ 2 + 2 1 ⋱ h ⎤ ⎥ ⎥ ⎥ ⎦ , = + -( ) 2 ( + ) -( + ) + -( ) 2 . ( 26
)
Finally, ̂ nl and ̂ ext represent respectively the Fourier coefficients of the vectors nl and ext . In the case of a harmonic excitation, the expression of ̂ ext is straightforward. However, it may be more difficult to determine the vector ̂ nl depending on the expression of the nonlinearity involved in the system. For geometrical nonlinearities (polynomial nonlinearities) such as those present in equation ( 16), an analytical expression can be derived for ̂ nl . Yet, this task may turn out to be cumbersome and the Alternating Frequency-Time (AFT) [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF] is generally preferred. This approach relies on expressing the vector ̂ in the time domain thanks to an inverse discrete Fourier transform, then evaluating the nonlinear forces nl ( ), after which the Fourier coefficients of the nonlinear forces in the frequency domain ̂ nl are derived using a direct Fourier transform. For harmonic analyses, the HBM scheme is generally coupled with a continuation method, to better capture some nonlinear effects such as turning points situations. In the context of the present work, a pseudo arc-length continuation method [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF] with a tangent prediction is implemented. The whole system including HBM and continuation equations is finally solved iteratively using a nonlinear solver [START_REF] Powell | A Fortran subroutine for solving systems of nonlinear algebraic equations[END_REF].

Branch-switching algorithm

When studying nonlinear geometric systems, abrupt qualitative changes in the solutions may appear [START_REF] Nayfeh | Applied nonlinear dynamics: analytical, computational, and experimental methods[END_REF]. Theses changes come from bifurcation phenomena such as turning or branching points. The former is generally handled well by the continuation procedure. The latter, however, cannot be overcome simply by classical continuation techniques since multiple branches of solutions appear simultaneously. Moreover, the challenge is all the greater as internal resonances often occur through new branches of solutions at the branching point bifurcation [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]. The branching analysis will be performed in three steps:

1. Detection of the branching point; 2. Localization of the branching point; 3. Calculation of the tangents to the branching solutions and tracking of the branches.

Since the branching point is by definition a singular point of the system (25), its detection is realized by monitoring the determinant of the Jacobian matrix of the system. A change of sign of the determinant indicates the presence of a bifurcation point.

The precise location of the branching point is computed by exploiting its singularity property. To this end, Seydel [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF] proposes to localize precisely the branching point by solving the following system for the unknown vector ̂ , , , :

( ̂ , , , ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ( ̂ , ) + , ̂ ( ̂ , ) , ( ̂ , ) -1 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = ( 27
)
where ( ̂ , ) refers to [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF], , ̂ and , are the derivatives of with respect to ̂ and , is a non-zero left eigenvector of the Jacobian matrix , ̂ , and a small scalar introduced to prevent singularity of the system. The last step of the branching analysis consists in calculating the tangents to the branching solutions at the branching point, noted * in the following. The tangents are obtained by means of the branch-switching algorithm introduced by Kuznetsov in [START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF]. They take the following form:

= 1 1 + 2 2 ( 28
)
with 1 and 2 two known vectors spanning the null space of the Jacobian matrix of system (25) and 1 and 2 two scalar solutions of the following quadratic form:

11 2 1 + 2 12 1 2 + 22 2 2 = 0 (29)
where the coefficients are defined by

= ⟨ * , ( i , j )⟩ = ∑ =1 +1 ∑ , =1 * 2 ( ) | | | | | = * ( i ) ( j ) , = 1, 2 (30)
In expression (30), = ̂ , , ( j ) is the m th component of the vector j , * represents a left eigenvector of the Jacobian matrix of evaluated at the bifurcation point * , and is the Hessian matrix of . Then, the calculation of the tangents is defined by the estimation of the couple 1 , 2 solution of the equation ( 29). An efficient procedure

225

for determining these couples is described in Appendix C.

Phenomenological insight model

In the present work, modal interactions between modes with commensurable frequencies will be examined. More specifically, attention will be paid to internal resonances between torsion and bending modes for which the frequency of the first mode is half that of the second mode. To facilitate the understanding of the mechanism driving the internal resonance, let us consider first the following nonlinear phenomenological system in modal form (2 modes only):

̈ 1 + 2 1 1 ̇̃ 1 + 2 1 ̃ 1 + 2 ∑ =1, ≥ ̃ ̃ + 2 ∑ =1, ≥ , ≥ ̃ ̃ ̃ = 0 ̈ 2 + 2 2 2 ̇̃ 2 + 2 2 ̃ 2 + 2 ∑ =1, ≥ ̃ ̃ + 2 ∑ =1, ≥ , ≥ ̃ ̃ ̃ = 0 cos ( ) (31)
with ̃ 1 and ̃ 2 the generalized coordinates, 1 and 2 the natural frequencies verifying the relation 2 ≃ 2 1 , 1 and 2 the damping ratios associated with the modes, and , , , the coefficients of the nonlinear quadratic and cubic terms. The excitation frequency is considered to be near 2 .

The system (31) is solved by means of the multiple scales method [START_REF] Nayfeh | Perturbation methods[END_REF]. The full development of the analysis is available in Appendix D. The procedure finally leads to two steady-state possibilities:

• the linear case,

1 = 0 and 2 = 2 2 2 2 + 2 2 1∕2 (32)
• the nonlinear case,

2 1 = Λ 1 ± 4 2 -Λ 2 2 1∕2 and 2 = 2 1 | | 12 | | 4 2 1 + 2 -1 2 1∕2 (33)
where 1 and 2 are directly linked to the modal amplitudes of ̃ 1 and ̃ 2 . The other terms in [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF] and [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF], as well as the expressions of Λ 1 and Λ 2 , refer to damping ( ), forcing ( ), and detuning ( ) terms and are introduced in Appendix D. The solutions [START_REF] Urabe | Galerkin's procedure for nonlinear periodic systems[END_REF] and ( 33) are very similar to those obtained by Nayfeh et al. in [START_REF] Nayfeh | Nonlinear coupling of pitch and roll modes in ship motions[END_REF]. For the linear solution (32), 1 = 0 implies that there is no internal resonance. However, the nonlinear solution (33) corresponds to a 1:2 internal resonance associated with the nonlinear terms 12 and 11 (see Appendix D). This observation directly implies that the appearance of a 1:2 internal resonance is only driven by these two nonlinear modal coefficients.

This example shows that even though all nonlinear coefficients may be involved in the nonlinear dynamic behavior, only a subset of them (here 12 and 11 ) have an impact on the existence of a 1:2 internal resonance. This remark must be kept in mind as it will be usefully referred to further in this paper.

Analysis and parametric studies

The previous sections have introduced a simplified blade model including pre-twist and pre-bending along with its discretized equation of motion. The different numerical methods to deal with the model have also been presented. The purpose of this section is to study the internal resonance involving the first torsional mode (noted 1T) and the second bending mode (noted 2B) of the structure. More specifically, the study aims to identify the physical nonlinear contributions responsible for the occurrence of internal resonance.

To this end, four different cases constructed from the blade model described previously are to be studied. These cases vary mostly according to their cross-section shape, pre-bending, and rotational speed. The underlying strategy is to isolate the effects of each contribution on the internal resonance. The parameters of the four configurations are given in Table 5. For each case, the geometrical properties are adjusted so that the frequency of the second bending mode is close to twice the frequency of the first torsional mode ( 2B ≃ 2 1T ). Note that the cross-section types are chosen in such a way that an analytical expression for the warping function is available (see Appendix A). This prerequisite insures a correct evaluation of the warping-dependent constants defined in (6).

In the following, the four studied cases are split into two groups: rotating and non-rotating ones. When the rotational speed is taken into account, commensurability between the modes is achieved by a particular rotational speed.

First case: rectangular cross-section without pre-bending

For this first case, the blade cross-section is assumed to be rectangular and there is no pre-bending. The pre-twist, however, is taken into account in such a way that the angle at the blade tip is 10°greater than the angle at the blade root. The frequency response of the blade is studied by imposing a localized harmonic force at the blade tip with an excitation frequency close to the frequency of the second bending mode (i.e close to twice the frequency of the first torsional mode). The HBM is used by setting h = 5 harmonics.

With the aim of observing the response of the first torsional mode near half the frequency of excitation, the fundamental frequency of the truncated Fourier series ( 23) is modified. The fundamental frequency is shifted so that the new fundamental frequency is equal to ∕2 ( = 2 × 2 ∕ the fundamental excitation period). This setup allows the first torsional mode to eventually respond through harmonic coefficients 1 and 1 .

A Rayleigh damping is introduced by setting the damping ratios 1T and 2B of the modes 1T and 2B both equal to 0.05 %. Unless explicitly specified, this choice of parameters is conserved for the whole set of simulations.

In the following, the displacement found with the HBM (by solving [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF] and then using the expansion ( 23)) is projected onto the linear modal basis of the system using the relation

̃ = with = n ( 34
)
in order to display directly the modal coordinates associated with the linear modes of interest.

The frequency responses of the modal coordinates ̃ 1T and ̃ 2B associated with the first torsional and second bending modes are represented in Fig. 6 for increasing forcing amplitudes. The frequency responses of the second bending mode (Fig. 6a) exhibit a softening behavior that increases with the excitation amplitude. This aspect is directly related to the presence of the nonlinearities in the model. A similar trend is observable for the frequency responses of the first torsional mode (Fig. 6b), yet with a modal amplitude about 100 times lower.

The relatively weak amplitude response of the first torsional mode results from a nonlinear coupling between the two modes. It has nothing to do, however, with an internal resonance. This is visible by referring to the time signals of the modal coordinates ̃ 1T ( ) and ̃ 2B ( ) plotted in Fig. 7 near the peak, depicted by [ ], of the frequency responses (Fig. 6). The oscillation frequency of the torsion mode (Fig. 7b) is indeed equal to twice that of the bending mode (Fig. 7a). The torsion mode is actually excited by a nonlinear harmonic of the bending mode and no internal resonance occurs in this situation.

The absence of internal resonance for this first blade model case is due to the nonlinear modal coefficients involved in the modal equations of the studied modes. Indeed, as mentioned in section 3.3, the internal resonance is essentially driven by the nonlinear modal coefficients: 12 and 11 . For this case, these modal coefficients can be found by projecting [START_REF] Hodges | Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[END_REF] on the linear modes using [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF] to obtain modal equations similar to [START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF] for both 1T and 2B modes ( 12 is hence the coefficient multiplying the term ̃ 1T ̃ 2B for the first equation and 11 is the coefficient linked to ̃ 2 1T for the second equation). Looking at the numerical values of the nonlinear coefficients present in the modal equations of modes 1T and 2F (illustrated in Fig. 8), it is straightforward that the two essential coefficients are negligible in comparison with some of the other coefficients.

An explanation for this observation could lie in the double symmetry of the blade cross-section. Indeed, this symmetry has a double impact. It generates uncoupled linear modes (see Fig. 9), but also implies a drastic simplification of

q u q v q w q φ Ψ 1T Ψ 2B -30
-20 -10 0 In view of the conclusions drawn, the next section will focus on a structure whose double-symmetry property no longer holds in order to evaluate the impact of such an assumption.

Second case: airfoil cross-section

This section focuses on a blade with an airfoil cross-section, for which an analytical expression of the warping function is known (see appendix A). The section possesses only one single axis of symmetry and the parameters of the model are listed in Table 5. Again, no pre-bending and rotation are accounted for in this model and the pre-twist is similar to the one used before. As for the previous section, the frequency response near the second bending mode is investigated. coupling leads to the emergence of two solution branches on both frequency responses (Figs. 10a and10b). These branches originate from branching points and have been obtained thanks to the algorithm presented in section 3.2. The first bifurcated branch appears with near twice the frequency of the first torsional mode ( ∕ 1T ≃ 2), while the second branch shows up in the neighborhood of the second bending mode frequency ( ∕ 1T ≃ 2.01). A similar behavior has been observed in [START_REF] Palma | Nonlinear harmonic analysis of a blade model subjected to large geometrical deflection and internal resonance[END_REF] and is directly related to a 1:2 internal resonance between the first torsional and the second bending modes. This is verified by looking at the time signals of the modal coordinates ̃ 1T ( ) and ̃ 2B ( ) for the two bifurcated branches (Fig. 11).

For each branch, the bending mode (Figs. 11a and11c) oscillates at a frequency equal to the excitation frequency. Conversely, the oscillation frequency of the torsion mode (Figs. 11b and11d) is equal to half the excitation frequency, which by construction is very close to the first torsional mode frequency. Mode 1T is therefore resonating when mode 2B is directly excited. This proves that a 1:2 internal resonance happens.

Designing structures that exhibit such internal resonances can jeopardize their integrity. Among the most dramatic features of the internal resonance, note that the amplitude of the first torsional mode can be larger than the amplitude of the second bending mode, as seen in the first bifurcated branch (Fig. 10b).

Aiming to understand the reasons for the appearance of internal resonance, the nonlinear modal coefficients involving the 1T and 2B modes are represented in Figure 12 (as was done previously) so as to evaluate their impact on the nonlinear behavior. The two coefficients 12 and 11 are much higher this time than in the previous case. More specifically, both 12 and 11 (respectively in front of ̃ 1T ̃ 2B and ̃ 2 1T ) are derived through differentiation from a coefficient, let us call it ̃ 112 , in front of ̃ 2

1T ̃ 2B in the potential energy [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF]. Thus, understanding the origins of 12 and 11 comes down to identifying the parameters that strongly contribute to the value of ̃ 112 . For this purpose, one can rewrite , , , (and their derivatives) in the potential energy [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF] in terms of ̃ 1T and ̃ 2B by using ( 13) and [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF], and identify the coefficient in front of ̃ 2 1T ̃ 2B . After identification, all the contributions due to ̃ 2 1T ̃ 2B are synthesized in Figure 13. The rows of the table correspond to the contribution to the term ̃ 2B , while the columns refer to the ̃ 2 1T contribution.

The black square in Figure 13 shows that ̃ 112 essentially stems from a term involving the monomial (including its derivatives) in the potential energy. After calculation, it appears that this contribution is related to the term 2 ′ ′′ ′′ cos(2 ) in eq. ( 5). Therefore, the internal resonance for the airfoil cross-section blade is mainly driven by the nonlinear part of the in-plane bending moment induced by the torsion component and the blade twist angle . Logically, this term involves both bending and torsion components.

0 0.2 0.4 0.6 0.8 1 -2 -1 0 1 2 •10 -2
Normalized time t/T Modal amplitude q2B

(a) 0 0.2 0.4 0.6 0.8 -2 -1 0 1 2 •10 -2
Normalized time t/T Modal amplitude q1T

(b) 0 0.2 0.4 0.6 0.8 1 -2 -1 0 1 2 •10 -2
Normalized time t/T Modal amplitude q2B

(c) 0 0.2 0.4 0.6 0.8 -2 -1 0 1 2 •10 -2
Normalized time t/T Modal amplitude q1T The modes also play an essential role in the appearance of internal resonance. In fact, the symmetry-breaking of the cross-section strongly changes the mode shapes of the structure and the decoupling observed for the rectangular cross-section no longer exists for the airfoil cross-section (see Fig. 14). This directly impacts the magnitude of the coefficient ̃ 112 as it is obtained through modal projection.

(d)
The symmetry-breaking of the cross-section has served to highlighted one factor responsible for the appearance of internal resonance. Other contributions, such as pre-bending, may lead to the same phenomenon implying the same modes. This case is discussed in the next section. In this section, the rectangular cross-section blade is considered once again, but a pre-bending is now added. Taking a curved geometry into account leads to consideration of additional linear and nonlinear terms in the model. Additionally, this new design of the blade makes it possible to activate other nonlinear terms which dissymetrize the structure.

q u q v q w q φ Ψ 1T Ψ 2B -5 0
The frequency responses of the modal coordinates ̃ 1T and ̃ 2B are obtained using the same techniques as before and are represented in Figure 15 for a forcing amplitude equal to 10 N. A significant difference exists between the frequency responses obtained without pre-bending (Fig. 6) and those obtained with pre-bending (Fig. 15). Among the main differences, one can observe that pre-bending changes the behavior of the model from softening to hardening.

Moreover, similarly to the previous case, one can distinguish the presence of a bifurcated branch in addition to the main solution branch. A detailed analysis on the associated time responses, not provided here, would give similar results to those observed in Figure 11 and would prove an occurrence of 1:2 internal resonance. Again, the nonlinear modal coefficients 12 and 11 reach higher levels than for the first blade case (Fig. 16) and explain the appearance of internal resonance. As explained and performed previously, the objective is now to identify the origins of the internal resonance for the curved model. To this end, a similar approach as that employed in section 4.2 gives the results shown in Figure 18, illustrating all the possible contributions to the coefficient ̃ 112 associated with the cubic term ̃ 2 1T ̃ 2B in the potential energy [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF].

α
Two major contributions now stand out. They correspond to the monomials 2 and (including their derivatives). Referring to [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF] and after identification, these monomials are related to the terms ′ -′′2 sin(2 ) and ′ 2 ′′ ′′ cos [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF] . Once again, the internal resonance is driven by nonlinear parts of the in-plane bending moment induced by the torsion . Moreover, its magnitude is still controlled by the twist angle and the material and geometrical ′ properties. As in the previous case, the important terms imply both bending and torsion components.

As we can observe, the pre-bending term is not involved in any of the two identified contributions. However, its impact on the magnitude of the coefficient ̃ 112 occurs for the modes of the structure, as one can see by comparing Fig. 9 and Fig. 17. Unlike in Fig. 9, the first torsional and second bending modes of Fig. 17 are not decoupled and share common components of displacements. As explained previously, this directly impacts the intensity of ̃ 112 , whose calculation is based on modal projection.

Finally, when comparing Figures 13 and18, one can observe that the highest contribution obtained with prebending is about twice as large as the highest contribution obtained with the airfoil cross-section. Hence, it seems that pre-bending has a greater influence than symmetry-breaking of the cross-section on the magnitude of ̃ 112 , and thus on the appearance of internal resonance.

q u q v q w q φ Ψ 1T Ψ 2B - 5
0

Fourth case: rotating blade

The previous sections studied the internal resonance of the blade model for various configurations. The physical contributions allowing the internal resonance to happen were also highlighted. Each model was tuned so that the frequency of the second bending mode was close to twice the frequency of the first bending mode in a non-rotating situation.

In this section, the objective is to investigate the case of an internal resonance occurring at a constant rotating speed, a feature that has never been addressed in the literature. For this purpose, the blade is tuned so that commensurability between the frequencies will be obtained for a known rotating speed. For the sake of simplicity, the blade is similar to the one studied in section 4.3 and is then pre-curved with a rectangular cross-section. The parameters used to build the model are listed in Table 5 and the associated Campbell diagram (Fig. 19) displays a commensurability condition for a rotating speed close to 2540 RPM. The frequency responses of the modal coordinates ̃ 2B and ̃ 1T for both Ω = 0 RPM and Ω = 2540 RPM are represented respectively in Figures 20 and21. For each case, the excitation frequency is close to the second bending mode frequency at the corresponding rotating speed and normalization is performed by considering 1T for Ω = 2540 RPM. At rest (Ω = 0 RPM), the response is almost exclusively composed of the second bending mode. The response of the first torsional mode is negligible in comparison. In this situation, no internal resonance is noticeable since the frequencies of 1T and 2B modes are not commensurate. However, the frequency responses obtained for the rotational speed [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF]. Additionally, bifurcation branches are visible again and are related to an internal resonance. The time signals on the branches are not given here but would confirm the internal resonance as in the previous sections. This last case has investigated the situation of an internal resonance for a nonzero rotating speed. In this case, the dynamics of the system involves both linear and nonlinear stiffnesses but also gyroscopic, spin-softening, and centrifugal pre-load effects. This example proves that internal resonance may also arise from nonlinearities in rotating condition, even if linear frequencies are not commensurate at zero rotating speed.

Conclusion

This paper has given some insights into the coupling between the torsional and bending modes of a blade due to the existence of an internal resonance. The model developed consisted of a rotating pre-bended and twisted blade whose equations of motion are obtained through a Rayleigh-Ritz discretization. A 1:2 internal resonance between the first torsion mode and the second bending modes has been investigated using the harmonic balance method coupled with a branch-switching algorithm. The forced responses of four parametrized configurations of the blade have been studied.

It has been shown that no internal resonance is visible when dealing with a twisted rectangular cross-section blade, because of its geometrical symmetries. Furthermore, internal resonances have been observed when breaking the double symmetry of the cross-section and considering an airfoil-like cross-section. In this case, we have demonstrated that the internal resonance is driven by a unique coupling term involving bending and torsion components of the first torsional and second bending modes. It has also been identified that pre-bending has a large impact on the appearance of internal resonance. This impact seems even greater than that observed for the airfoil cross-section. Finally, a rotating blade has been investigated involving gyroscopic, spin-softening, and centrifugal pre-load effects. In this last example, an internal resonance has been also observed at the operating point corresponding to 2540 RPM.

The results exhibited in this paper can be used as general guidelines for design. They give some insights into how to avoid or control the occurrence of internal resonance and can help engineers in their strategies to achieve their performance objectives. Note that the developed model can also come in handy to study further nonlinear phenomena such as instabilities, multi-mode responses, etc. This appendix furnishes the analytical expression of the transformation matrix relating global to local axes. Based on [START_REF] Peters | Technical notes: The effects of second order blade bending on the angle of attack of hingeless rotor blades[END_REF] and considering the pre-bending , its expression is

= ⎡ ⎢ ⎢ ⎢ ⎣ 1 + + + + -

C. Useful method for calculating the tangents in the branch-switching algorithm

This appendix proposes a method to derive the tangents of the solution branches at a bifurcation point. As seen in section 3.2, the method comes down to finding the solutions of the quadratic form [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. Hence, here we suggest a procedure to solve that equation.

Equation ( 29) is first rewritten as the following matrix equation: (50)

Introducing (50) in (49) and setting = , one can retrieve the following relation:

1 = | | 2 | | (51)
where = ± √ 1 ∕ 2 . At this point, one equation for the two new unknowns 1 and 2 is obtained. Prescribing the norm of the sought tangents equal to Δ yields a second equation to the system. Different solutions then exist according to the signs of 1 and 2 . These are enumerated in Table 6.

In practice, cases 2 and 3 in Table 6 give the same tangent with the opposite direction. The same remark holds for cases 1 and 4. The correct direction to follow can be found using a scalar product between the computed tangent and the tangent evaluated at the previous iteration of the continuation procedure.

Once 1 and 2 are computed, 1 and 2 are deduced using = . Finally, the tangents are obtained using the relation [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF].

D. Details of the multiple scales analysis

This appendix details the developments of the multiple scales analysis used in section 3.3. For that purpose, the different time scales = are introduced, with a small dimensionless parameter. Taking into account the assumptions made in section 3.3, we define 2 1 = 2 + 1 and = 2 + 2 (52) with 1 and 2 two detuning parameters. The following orders of magnitude are also set:

0 = 2 , 1 1 = 1 and 2 2 = 2 (53)
Following the multiple scales method, expansions of ̃ 1 and ̃ 2 are of the form ̃ 1 = ̃ 11 ( 0 , 1 ) + 2 ̃ 12 ( 0 , 1 ) ̃ 2 = ̃ 21 ( 0 , 1 ) + 2 ̃ 22 ( 0 , 1 )

By expressing the time derivatives with respect to the new time scales and introducing the expansions (54) in equations [START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF]

Figure 1 :

 1 Figure 1: Pre-bent and twisted blade model in its initial and deformed states

Figure 2 :

 2 Figure 2: Blade cross-section at blade tip before and after deformation

Figure 3 :

 3 Figure 3: Campbell Diagram of the second validation step structure: proposed blade model [ ], Ansys [ ]

Figure 4 :

 4 Figure 4: Mode shapes of the second bending (2B) and first torsional (1T) modes obtained for Ω = 0 RPM

Figure 5 :

 5 Figure 5: Mode shapes of the second bending (2B) and first torsional (1T) modes obtained for Ω = 5000 RPM

Figure 6 :

 6 Figure 6: Frequency responses with F= [50 N, 100 N, 150 N, 200 N] of the second bending mode (a) and the first torsional mode (b) for the first blade model case

Figure 7 : 6 α

 76 Figure 7: Time signals (a): ̃ 2B , and (b): ̃ 1T at the resonant peaks [ ] of Figure 6

Figure 8 :

 8 Figure 8: Nonlinear modal coefficients associated with modes 1T and 2F in the two modal equations describing the first blade model case

Figure 9 :

 9 Figure 9: Eigenvector amplitudes in log scale of the first torsional mode 1T and the second bending mode 2B of the first blade model case.

  295

N

  . Di Palma et al.: Preprint submitted to Elsevier

Figure 10 Figure 10 :

 1010 Figure10represents the frequency responses of the first torsional mode and the second bending mode of this testcase. The two observed modes respond to the excitation near 2B , which shows the presence of mode coupling. This

Figure 11 :

 11 Figure 11: Time signals ̃ 2B and ̃ 1T at the resonant peaks [ ] of the first bifurcated branch (a) and (b), and of the second bifurcated branch (c) and (d) of figure 10

α 11 β

 11 11 α 12 β 12 α 22 β 22 α 111 β 111 α 112 β 112 α 122 β 122 α 222 β 222

Figure 12 :

 12 Figure 12: Nonlinear modal coefficients associated with modes 1T and 2F in the two modal equations describing the second blade model case

4 Figure 13 :

 413 Figure 13: Origins of the coefficient ̃ 112 associated with the cubic term ̃ 2 1T ̃ 2B in the potential energy of the second blade model case

Figure 14 : 4 . 3 .

 1443 Figure 14: Eigenvector amplitudes in log scale of the first torsional mode and the second bending mode of the second blade model case

Figure 15 :

 15 Figure 15: Frequency responses of the second bending mode (a) and the first torsional mode (b) for the third blade model case.

Figure 16 :

 16 Figure 16: Nonlinear modal coefficients associated with modes 1T and 2F in the two modal equations describing the third blade model case

Figure 17 :

 17 Figure 17: Eigenvector amplitudes in log scale of the first torsional mode and the second bending mode of the third blade model case

Figure 18 :

 18 Figure 18: Origins of the coefficient ̃ 112 associated with the cubic term ̃ 2 1T ̃ 2B in the potential energy of the third blade model case

Figure 19 :

 19 Figure 19: Campbell diagram of the fourth test case, highlighting the commensurability between the modes 1T and 2B near 2540 RPM

Figure 20 :

 20 Figure 20: Frequency responses of the second bending mode (a) and the first torsional mode (b) for the fourth blade model case (Ω = 0 RPM).

Figure 21 :

 21 Figure 21: Frequency responses of the second bending mode (a) and the first torsional mode (b) for the fourth blade model case (Ω = 2540 RPM).

  the matrix is invertible and there are an orthogonal matrix and a diagonal matrix such that

Table 3

 3 Comparison of the natural frequencies of the straight beam

	Modes	Euler-Bernoulli theory [Hz]	Model [Hz]	Relative error [%]
	1	30.71	30.71	8.13 × 10 -9
	2	185.25	185.25	1.87 × 10 -9
	3	192.46	192.46	1.39 × 10 -3
	4	307.11	307.11	8.13 × 10 -9
	5	538.89	539.06	3.09 × 10 -2
	6	555.75	555.75	6.87 × 10 -4
	7	926.25	927.58	1.44 × 10 -1

Table 4

 4 Parameters of the twisted and curved model for the second validation step

		Geometry			Material			Twist		Pre-bending
	L [m]	w [m]	h [m]	[kg m -3 ]	E [GPa]			( )		( )
	1.3	0.4	0.04	1550	100	0.3	2	+	6	0.1 2

Table 5

 5 Parameters of the four blade model cases

	Cases	Geometry				Material				Twist	Pre-bending Rotational speed
		Section type L [m] w [m] h [m]	[kg m -3 ] E [GPa]			( )		( )	Ω [RPM]
	1	rectangle	1.3	0.97 0.065	1550	3.8	0.4	2	+	18	none	0
	2	airfoil	1.2	1.1	0.12	1550	1.6	0.4	2	+	18	none	0
	3	rectangle	1.2	0.96 0.065	1550	3.8	0.4	2	+	18	0.1 2	0
	4	rectangle	1.2	0.6	0.06	1550	10	0.4	2	+	18	0.1 2	2550

  11 β 11 α 12 β 12 α 22 β 22 α 111 β 111 α 112 β 112 α 122 β 122 α 222 β 222

	Amplitude [SI]	10 2 10 4 10 6
		Nonlinear modal coefficients

  In the absence of deformation and pre-bending, one can verify that reduces to a rotation matrix of angle about the global axis x.Using the ordering scheme introduced in section 2.2 and assuming a moderate pre-bending, it reads Hence, to the second order, no distinction is made between the curvilinear coordinate and the coordinate . This remark implies that d.∕ d ≃ d.∕ d . Using this simplification, the transformation matrix expanded to second order becomes

	d d	=	1 + √	1	d d	2	≃ 1	(46)
			⎡ ⎢			1 -	1 2		′ 2 -	1 2	′ + ′ 2	′	′ + ′	⎤ ⎥
		=	⎢ ⎢ ⎢	-′ cos( + ) + ( ′ + ′ ) sin( + ) cos + ̂ + ′ ( ′ + ′ ) 1 -	1 2	′ 2	sin( + ̂ ) 1 -	1 2	′ + ′ 2	⎥ ⎥ ⎥	(47)
			⎢ ⎢ ⎣	[ ′ sin( + ) -( ′ + ′ ) cos( + )] -sin( + ̂ + ′ ( ′ + ′ )) 1 -	1 2	′ 2	cos( + ̂ ) 1 -	1 2	( ′ + ′ )	2	⎥ ⎥ ⎦
	noting (.) ′ = d.∕d and with,
	̂ ≡ -∫ 0		′′ ( ′ + ′ ) d	(48)
				1 2	1 2	1 sin ̄ + + + + cos ̄ + -1 cos ̄ + + + + sin ̄ + -1 1 2 2	1 cos ̄ -sin ̄ + + + + 1 sin ̄ + cos ̄ + + + +	2 sin ̄ 2 cos ̄	⎤ ⎥ ⎦ ⎥ ⎥	(44)
	where,
	1 =	√ 1 -+ 2 -+ + + 2
	2 = ̄ = + -∫ 0 √ 1 -+ 2	++ 1 -+ 2 -+ + + 2 + + + √	d -∫ 0	1 -+ + + 2 + + + 2 ++ + ++ + 1 -+ 2 -+ + + 2 √	d	(45)

Table 6

 6 CaseAssumption

			1		2
	1	> 0 and 2 > 0	2		Δ
	2	> 0 and 2 < 0 -2		‖ ‖ ‖ 11 + 12 1 + -Δ	21 + 22 2	‖ ‖ ‖
	3	< 0 and 2 > 0	2	‖ ‖ ‖	-11 + 12 1 + -21 + 22 2 Δ	‖ ‖ ‖
	4	< 0 and 2 < 0 -2		‖ ‖ ‖ 11 + 12 1 + -Δ	21 + 22 2	‖ ‖ ‖
				‖ ‖ ‖	-11 + 12 1 + -21 + 22 2	‖ ‖ ‖

  ̃ 12 = -2 0 1 ̃ 11 -2 1 0 ̃ 11 + 11 ̃ 11 ̃ 11 + 12 ̃ 11 ̃ 21 + 22 ̃ 21 ̃ 21 2 0 ̃ 22 + 2 2 ̃ 22 = -2 0 1 ̃ 21 -2 1 0 ̃ 21 + 11 ̃ 11 ̃ 11 + 12 ̃ 11 ̃ 21 + 22 ̃ 21 ̃ 21 + cos ( )

		, we get
	• for the first order
			2 0 ̃ 11 + 2 1 ̃ 11 = 0 2 0 ̃ 21 + 2 2 ̃ 21 = 0	(55)
	• for the second order 2
			2 0 ̃ 12 + 2 1 (56)
	where	= ∕ .
	The solutions of equations (55) have the following form:
		̃ 11 = 1 ( 1 ) i 1 0 + . ̃ 21 = 2 ( 1 ) i 2 0 + .	(57)
	where . stands for complex conjugate. Injecting solutions (57) into the equations (56) yields the following solvability
	conditions:
	⎧ ⎪ ⎨ ⎪ ⎩	-2 1 -2 2	′ 1 + 1 1 -12 2 ̄ 1 ′ 2 + 2 2 -11 2 1 i 1 1 + -i 1 1 = 0 2 i 2 1 = 0	(58)

Acknowledgements

The authors would like to acknowledge the financial support of Safran Aircraft Engines.

Ω = 2540 RPM are drastically different. Notice the softening behavior of the second bending mode when increasing

Appendix

A. Analytical form of the warping function for rectangular and airfoil cross-sections

This appendix gives a brief reminder of the warping function calculation and details the exact analytical expressions of the warping function for rectangular and airfoil cross-sections. Most of the developments presented here could be found in greater detail in [START_REF] Timoshenko | Theory of Elasticity[END_REF].

The warping function of a prismatic bar of section under torsion is the solution of the following system:

For a rectangular cross-section of sides 2 × 2 , one can verify that the following expression is the solution of problem [START_REF] Powell | A Fortran subroutine for solving systems of nonlinear algebraic equations[END_REF] [START_REF] Timoshenko | Theory of Elasticity[END_REF] ( ,

Equation ( 36) represents the analytical exact warping function of the rectangular cross-section. For practical purposes the infinite sum can be truncated, keeping only 5 terms since the convergence is very fast.

For more complex cross-sections, solving system (35) may be impractical. To simplify the procedure, a useful alternative consists in solving the following equivalent system:

where Φ is the Prandtl stress function related to the warping function by the following expressions:

In some applications, it is easier to solve system (37) in order to derive the expression of the warping function. Saint Venant showed that the system (37) admits an infinite number of solutions of the general form [START_REF] De Saint-Venant | Mémoire sur la torsion des prismes: avec des considérations sur leur flexion ainsi que sur l'équilibre intérieur des solides élastiques en général: et des formules pratiques pour le calcul de leur résistance à divers efforts s'exerçant simultanément[END_REF],

Based on these general solutions (39), Wang [START_REF] Wang | Exact torsion solutions for certain airfoil shapes[END_REF] proposes to choose the three-term sum in even power whos expression is derived as follows:

By adjusting the constants 0 , 2 and 4 , the application of the second equation of system [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF] can lead directly to the boundary shape equation of an airfoil cross-section of chord :

with

The parameter 0 in eq [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF] has been added in this paper to shift the airfoil cross-section to force the shear center to be located at the intersection of the principal axes of the cross-section. Using [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF], the expression of the warping function for the airfoil cross section is straightforward and becomes

B. Expression of the transformation matrix T

where primes denote differentiation with respect to 1 . One then uses the notation

and

Introducing ( 59) and (60) in the solvability conditions (58) and splitting real and imaginary parts give (62)

• the nonlinear case,