Marielle Debeurre
email: marielle.debeurre@ensam.eu

Aurélien Grolet
email: aurelien.grolet@ensam.eu

Bruno Cochelin
email: bruno.cochelin@centrale-marseille.fr

Olivier Thomas
email: olivier.thomas@ensam.eu

Finite element computation of nonlinear modes and frequency response of geometrically exact beam structures

An original method for the simulation of the dynamics of highly flexible slender structures is presented. The flexible structures are modeled via a finite element (FE) discretization of a geometrically exact two-dimensional beam model, which entirely preserves the geometrical nonlinearities inherent in such systems where the rotation of the cross-section can be extreme. The FE equation is solved by a combination of harmonic balance (HBM) and asymptotic numerical (ANM) methods. The novel solving scheme is rooted entirely in the frequency domain and is capable of computing both the structure's frequency response under periodic external forces as well as its nonlinear modes. An overview of the proposed numerical strategy is outlined and simulations are shown and discussed in detail for several test cases.

Introduction

This paper considers the numerical simulation of the dynamical behaviors of highly flexible slender beam structures subjected to geometrical nonlinearities. The geometrical nonlinearities are observed here due to the thinness of the structure, which leads to a very small stiffness in the transverse direction of deformation. Large displacements are therefore possible even for moderate external forces. This is especially true of structures with boundary conditions analogous to a cantilever (Fig. 1(b)), for which the main source of geometrical nonlinearities is not a membrane/bending coupling -encountered for instance in clamped-clamped beams, plates and shells (Fig. 1(a)) -but comes from the rotation of the cross-section, which can be very large (see e.g. [START_REF] Touzé | Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques[END_REF][START_REF] Lacarbonara | Refined models of elastic beams undergoing large in-plane motions: Theory and experiment[END_REF]). This latter case is associated with moderate nonlinear effects, which become significant for very large displacements of the structure (on the order of its characteristic length). Proper simulation of such extreme deformations of various beam structures is precisely the purpose of this article.

This field is rich with industrial engineering applications (see the introduction of [START_REF] Paï | Highly Flexible Structures: Modeling, Computation, and Experimentation[END_REF]), for example where large amounts of wiring or cable-bundling are necessary, such as in the automotive industry. Recent work has addressed this topic in statics [START_REF] Cottanceau | A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables[END_REF] and dynamics [START_REF] Linn | Kinetic aspects of discrete cosserat rods based on the difference geometry of framed curves[END_REF][START_REF] Lang | Multi-body dynamics simulation of geometrically exact Cosserat rods[END_REF][START_REF] Brüls | Lie group generalized-𝛼 time integration of constrained flexible multibody systems[END_REF][START_REF] Sonneville | Geometrically exact beam finite element formulated on the special Euclidean group SE(3)[END_REF]. Other applications of nonlinear beam models are the design of rotating blades, such as in helicopters, turbomachinery or wind turbines [START_REF] Hodges | Nonlinear Composite Beam Theory[END_REF][START_REF] Hodges | Geometrically-exact, intrinsic theory for dynamics of curved and twisted anisotropic beams[END_REF][START_REF] Palacios | Nonlinear normal modes in an intrinsic theory of anisotropic beams[END_REF][START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF], highly flexible structures in aerospace such as high-altitude, long endurance aircraft [START_REF] Su | Strain-based geometrically nonlinear beam formulation for modeling very flexible aircraft[END_REF][START_REF] Perroni | Effects of geometric nonlinearity on flexible wing structures[END_REF], vibrating wings of nano-drones [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF][START_REF] Faux | Dynamic simulation and optimization of artificial insect-sized flapping wings for a bioinspired kinematics using a two resonant vibration modes combination[END_REF], flexible cables for cable car systems [START_REF] Gerstmayr | A 3D finite element method for flexible multibody systems[END_REF][START_REF] Gerstmayr | A geometrically exact beam element based on the absolute nodal coordinate formulation[END_REF][START_REF] Bertrand | A robust and efficient numerical finite element method for cables[END_REF], vibration isolation with pinched loops [START_REF] Virgin | Vibration isolation using extreme geometric nonlinearity[END_REF] and micro/nanosystems such as mass sensors and atomic force microscopes [START_REF] Dick | Parametric resonance and pattern selection in an array of microcantilevers interacting through fringing electrostatic fields[END_REF][START_REF] Rabenimanana | Mass sensor using mode localization in two weakly coupled MEMS cantilevers with different lengths: Design and experimental model validation[END_REF][START_REF] Ben Shaya | Automatic locking of a parametrically resonating, base-excited, nonlinear beam[END_REF][START_REF] Vincent | Nonlinear polarization coupling in freestanding nanowire/nanotube resonators[END_REF]. Another area of interest in this vein of research is the exploitation of nonlinear cantilever beams for piezoelectric control and energy-harvesting applications, see e.g. [START_REF] Elahi | Piezoelectric Aeroelastic Energy Harvesting[END_REF][START_REF] Mahmoodi | Non-linear vibrations and frequency response analysis of piezoelectrically driven microcantilevers[END_REF][START_REF] Mahmoodi | Subharmonics analysis of nonlinear flexural vibrations of piezoelectrically actuated microcantilevers[END_REF][START_REF] Yang | On the efficiency of piezoelectric energy harvesters[END_REF][START_REF] Leadenham | Unified nonlinear electroelastic dynamics of a bimorph piezoelectric cantilever for energy harvesting, sensing, and actuation[END_REF][START_REF] Guillot | Theoretical and experimental investigation of a 1:3 internal resonance in a beam with piezoelectric patches[END_REF].

A common model for numerical computation related to beam structures is the so-called ''inextensible beam'' model of Crespo da Silva and Glynn [START_REF] Crespo Da Silva | Nonlinear flexural-flexural-torsional dynamics of inextensible beams. I. equations of motion[END_REF][START_REF] Crespo Da Silva | Nonlinear flexural-flexural-torsional dynamics of inextensible beams. II. forced motions[END_REF][START_REF] Nayfeh | Linear and Nonlinear Structural Mechanics[END_REF]. In this model the length of the beam is supposed to be a constant, which allows the axial motion to be condensed into the transverse motion, thereby reducing the dynamics to a single degree of freedom partial differential equation. In addition, the geometrical nonlinearities are truncated up to the third order, retaining only cubic nonlinearities [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF]. This model served as the basis for a huge number studies, but is limited to moderately large amplitudes due to the truncation of the nonlinear terms.

Highly flexible structures, however, can be subjected to deflections of extreme amplitude that necessitate richer beam models in order to be simulated appropriately. One solution is to turn to classical three dimensional finite element simulations [START_REF] Zienkiewicz | The Finite Element Method[END_REF], available in most commercial FE codes, which incorporate geometrical nonlinearities and do not impose restrictions on the maximum deflection that can be simulated accurately [START_REF] Touzé | Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques[END_REF][START_REF] Shen | Comparison of reduction methods for finite element geometrically nonlinear beam structures[END_REF][START_REF] Vizzaccaro | High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point[END_REF][START_REF] Farokhi | Three-dimensional nonlinear extreme vibrations of cantilevers based on a geometrically exact model[END_REF]. But, in this case, an enormous number of elements (and, consequently, of degrees of freedom, leading to sometimes prohibitive computation times) is required for accurate simulation of the slender structures considered in this article, largely due to the small thickness and their flexibility in bending. This justified the development of an alternative method, the so-called geometrically exact beam model, introduced by Reissner [START_REF] Reissner | On one-dimensional finite-strain beam theory: The plane problem[END_REF][START_REF] Reissner | On finite deformations of space-curved beams[END_REF] and Simo [START_REF] Simo | A finite strain beam formulation. The three-dimensional dynamic problem. Part I[END_REF]. This model considers beam kinematics without truncation of the nonlinear terms and, thus, does not restrict the amplitude of the rotation of the cross-section [START_REF] Meier | Geometrically exact finite element formulations for slender beams: Kirchhoff-Love theory versus Simo-Reissner theory[END_REF]. The model is therefore able to capture nonlinear behavior accurately even at extreme displacements. To reduce the computation time, reduced order models can also be used, see e.g. the recent review work of Touzé et al. [START_REF] Touzé | Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques[END_REF] and some applications [START_REF] Shen | Comparison of reduction methods for finite element geometrically nonlinear beam structures[END_REF][START_REF] Vizzaccaro | High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point[END_REF][START_REF] Lülf | An integrated method for the transient solution of reduced order models of geometrically nonlinear structures[END_REF].

A known difficulty in putting geometrically exact beam models into practice is linked to the treatment of the rotations for threedimensional -axial/bending/torsion -deflections of the beam. Different methods for parameterization of the rotations (rotation matrices, quaternions. . . ) have been proposed by Simo and Vu-Quoc [START_REF] Simo | A three-dimensional finite-strain rod model. Part II: Computational aspects[END_REF], Cardona and Géradin [START_REF] Cardona | A beam finite element non-linear theory with finite rotations[END_REF], Ibrahimbegović [START_REF] Ibrahimbegović | On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements[END_REF], Crisfield and Jelenić [START_REF] Crisfield | Objectivity of strain measures in geometrically exact 3D beam theory and its finite element implementation[END_REF][START_REF] Jelenić | Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics[END_REF] and Zupan et al. [START_REF] Zupan | The quaternion-based three-dimensional beam theory[END_REF], among others. These studies aim to implement a numerical finite element-based discretization of the model, allowing for the simulation of structures with arbitrarily complex geometries obtained through the assembly of beam elements.

In parallel, a series of recent works targets the dynamics of a single cantilever beam in the frequency domain, using the same geometrically exact model with an inextensibility constraint but without any truncation of the geometrical nonlinearities. This method leads to a simple formulation in the form of partial differential equations with transverse degrees of freedom only, which are discretized by an expansion onto a normal mode basis. The strategy produces highly accurate results when compared with reference numerical FE simulations or experiments [START_REF] Farokhi | Three-dimensional nonlinear extreme vibrations of cantilevers based on a geometrically exact model[END_REF][START_REF] Farohki | Geometrically exact extreme vibrations of cantilevers[END_REF][START_REF] Farohki | Experimentally validated geometrically exact model for extreme nonlinear motions of cantilevers[END_REF][START_REF] Ghayesh | Extremely large dynamics of axially excited cantilevers[END_REF], although it is restricted to the case of a single cantilever beam.

Most numerical strategies able to solve the geometrically exact beam model, available either in commercial finite-element codes or addressed in the literature, are related to time integration due to their highly nonlinear nature. However, when studying the nonlinear dynamics of such highly flexible structures, a lot of complex vibratory phenomena may be observed, including bifurcations, instabilities, bistability, qualitative changes in the nature of the solutions, quasi-periodic and chaotic responses, energy transfers between vibration modes, etc. [START_REF] Vincent | Nonlinear polarization coupling in freestanding nanowire/nanotube resonators[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF][START_REF] Mangussi | Resonance in a vibrating beam: a zoo of nonlinear resonance peaks[END_REF][START_REF] Cusumano | Chaotic non-planar vibrations of the thin elastica, part 1: experimental observation of planar instability[END_REF]. In this case, simulations in the frequency domain are advantageous since they directly target the periodic steady-state vibrations under harmonic forcing and can be easily integrated into bifurcation and stability analysis tools. Moreover, the nonlinear resonance pattern of a given structure can be efficiently mapped by the so-called nonlinear normal modes [START_REF] Kerschen | Nonlinear normal modes, Part I: A useful framework for the structural dynamicist[END_REF], equivalently defined as families of periodic solutions of the conservative system or invariant manifolds of the phase space [START_REF] Touzé | Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques[END_REF][START_REF] Touzé | Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes[END_REF][START_REF] Shaw | Nonlinear normal modes and invariant manifolds[END_REF]. The nonlinear normal modes are often depicted as backbone curves in an amplitude/frequency (also known as a frequency response) plot. This paper proposes a novel numerical strategy to compute the nonlinear dynamics of highly flexible structures in the frequency domain. It is based on a finite element discretization of the geometrically exact beam model and computes periodic solutions using a combination of the harmonic balance method (HBM) and a numerical continuation tool based on the asymptotic numerical method (ANM). The numerical strategy is capable of calculating both the frequency response under harmonic forcing and the nonlinear normal modes of structures of arbitrary shape, as long as they can be discretized into beam elements. Due to the implementation of the continuation method for solving, both steady-state stable and unstable periodic solutions can be computed. In this article, only plane structures with in-plane deformations are considered; the extension to a full three-dimensional case is postponed to a future article. The methodology is presented in detail in Section 2, including a discussion on the numerical properties and efficiency of the strategy, with several test cases gathered in Section 3.

Geometrically exact beam model and solution

In this section, the fundamental equations governing the motion of flexible beam structures in two dimensions (2D) are presented. Previous work derives in detail the kinematics and the equations of motion corresponding to a geometrically exact beam model [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF]. For this reason, only a summary of the equations central to the argument of this article is presented. The interested reader may refer to Sections 2 and 4 of [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF] for a more detailed derivation.

Kinematics

A depiction of the system under investigation is shown in Fig. 2. As mentioned, the geometrically exact beam model is here restricted to two dimensions, such that the centerline of any beam element lies in the (𝐞 𝑋 , 𝐞 𝑌 ) plane of a global basis (𝐞 𝑋 , 𝐞 𝑌 , 𝐞 𝑍 ) and undergoes deformations in the same plane.

Based on Timoshenko kinematics for the beam, we assume that any cross-section of the beam undergoes a rigid-body motion. We consider a differential beam element of length d𝑥 with its centerline aligned with the 𝐞 𝑥 direction, where (𝐞 𝑥 , 𝐞 𝑦 , 𝐞 𝑍 ) is a local frame that may depend (in the case of a curved structure) on 𝑥, the curvilinear coordinate that defines the location of the cross-sections. The angle 𝜓(𝑥) = (ê 𝑋 , 𝐞 𝑥 ) defines its orientation. As a consequence, the displacement of any point 𝑃 in the cross-section is related to the displacement field of the centerline by the following relation:

𝐮(𝑥, 𝑦, 𝑡) = [𝑢(𝑥, 𝑡) -𝑦 sin 𝜃(𝑥, 𝑡)] 𝐞 𝑥 + [𝑤(𝑥, 𝑡) + 𝑦 (cos 𝜃(𝑥, 𝑡) -1)] 𝐞 𝑦 , (1)
where (𝑥, 𝑦) are the coordinates of 𝑃 in the frame (𝐞 𝑥 , 𝐞 𝑦 ), 𝑡 is the time, 𝑢(𝑥, 𝑡) and 𝑤(𝑥, 𝑡) are the axial and transverse displacements, respectively, of the corresponding point of the centerline 𝐶 at axial position 𝑥 and time 𝑡 and 𝜃(𝑥, 𝑡) is the rotation of the cross-section, as shown in Fig. 2.

The expression for the strains, incorporating large displacements, is given by the Green-Lagrange strain tensor, denoted here Ẽ, and expressed as (see e.g. [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF]):

Ẽ = 1 2 ( ∇𝐮 + ∇ T 𝐮 + ∇ T 𝐮 ∇𝐮 ) . ( 2
)
For the sake of efficiency, the exact expressions of the Green-Lagrange strain tensor as a function of 𝑢, 𝑤 and 𝜃 are not recopied here. Instead, a consistent linearization to the Green-Lagrange strains is computed in order to simplify the expressions without loss of accuracy. In this process, only the part of the strain terms relating purely to the local strains is linearized, maintaining the rigid-body rotation part fully exact. It is relevant so long as the local strains remain relatively small, which is the case here, in spite of large rigid-body displacements (the latter will be assumed throughout the remainder of the paper). The resulting simplified strain tensor, denoted 𝐄, is then given in the basis (𝐞 𝑥 , 𝐞 𝑦 , 𝐞 𝑍 ), by the following [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF][START_REF] Géradin | Flexible Multibody Dynamics: A Finite Element Approach[END_REF][START_REF] Felippa | Nonlinear Finite Element Methods, Chapter 9: The TL Timoshenko Plane Beam Element[END_REF]:

𝐄 = ⎛ ⎜ ⎜ ⎝ 𝑒 -𝑦𝜅 𝛾∕2 0 𝛾∕2 0 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , ( 3
)
where the quantities 𝑒, 𝛾 and 𝜅 represent, respectively, the axial strain, the shear strain and the beam curvature, defined as:

𝑒 = (1 + 𝑢 ′ ) cos 𝜃 + 𝑤 ′ sin 𝜃 -1, 𝛾 = 𝑤 ′ cos 𝜃 -(1 + 𝑢 ′ ) sin 𝜃, ( 4
)
𝜅 = 𝜃 ′ ,
where

• ′ = 𝜕•∕𝜕𝑥.
The stress tensor (energetically conjugated to the Green-Lagrange strain tensor) is given by the second Piola-Kirchhoff stress tensor, denoted 𝐒. In what follows, we consider that the material is homogeneous and isotropic and that a linear Kirchhoff-Saint-Venant constitutive relation is assumed between 𝐄 and 𝐒, so that the following relation holds:

𝑁 = ∫ 𝐴 𝑆 𝑥𝑥 d𝐴 = 𝐸𝐴𝑒, 𝑇 = ∫ 𝐴 𝑆 𝑥𝑦 d𝐴 = 𝑘𝐺𝐴𝛾, ( 5
)
𝑀 = ∫ 𝐴 𝑦𝑆 𝑥𝑥 d𝐴 = 𝐸𝐼𝜅,
where 𝑁, 𝑇 and 𝑀 are the axial force, shear force and bending moment, respectively, 𝐸 is the Young's modulus, 𝐺 the shear modulus, 𝐴 and 𝐼 the area and second moment of area of the cross-section, respectively, and 𝑘 the shear correction factor [START_REF] Hutchinson | Shear coefficients for Timoshenko beam theory[END_REF]. The generalized forces are grouped into a vector, denoted Ŝ:

Ŝ = [ 𝑁 𝑇 𝑀 ] T , ( 6
)
where T defines the transpose of the vector or matrix.

The equations of motion are derived in the weak form using the principle of virtual work, which is stated in [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF][START_REF] Lacarbonara | Nonlinear Structural Mechanics: Theory, Dynamical Phenomena and Modeling[END_REF] as, for all virtual displacement 𝛿𝐮:

𝛿𝑊 𝑎 + 𝛿𝑊 𝑖 = 𝛿𝑊 𝑒 , (7)
with the inertial (𝛿𝑊 𝑎 ), internal (𝛿𝑊 𝑖 ) and external (𝛿𝑊 𝑒 ) virtual work expressed as, respectively:

𝛿𝑊 𝑎 = ∫ 𝑉 𝜌 ü ⋅ 𝛿𝐮 d𝑉 , ( 8
)
𝛿𝑊 𝑖 = ∫ 𝑉 𝐒 ∶ 𝛿𝐄 d𝑉 , ( 9
)
𝛿𝑊 𝑒 = ∫ 𝜕𝑉 𝐟 ⋅ 𝛿𝐮 d𝐴 + ∫ 𝑉 𝐛 ⋅ 𝛿𝐮 d𝑉 , ( 10
)
where 𝑉 represents the domain of the beam in the reference configuration, 𝜕𝑉 its frontier, 𝜌 the density of the material, ü = 𝜕 2 𝐮∕𝜕𝑡 2 the acceleration of point 𝑃 , 𝛿𝐮 the virtual displacement vector and 𝛿𝐄 the first variation of the strain induced by 𝛿𝐮. In addition, 𝐟 and 𝐛 represent, respectively, the external surface and body forces.

Finite element model

In this section we present the discretization used to derive the finite element equation of motion. First, the beam is discretized into 𝑁 e elements of length 𝐿 𝑒 , that may depend on 𝑥 in the case of a non-uniform mesh. The element used here is the classical two-node Timoshenko beam element with linear shape functions [START_REF] Felippa | Nonlinear Finite Element Methods, Chapter 9: The TL Timoshenko Plane Beam Element[END_REF]. At each node 𝑖 of an element, there are three degrees of freedom (𝑢 𝑖 , 𝑤 𝑖 , 𝜃 𝑖 ) corresponding to the axial displacement, transverse displacement and rotation of the cross-section, respectively. For a given element, we introduce the vector 𝐪 𝑒 containing the degrees of freedom at the two nodes of the element (denoted node 1 and node 2):

𝐪 𝑒 = [ 𝑢 1 𝑤 1 𝜃 1 𝑢 2 𝑤 2 𝜃 2 ] T . ( 11
)
We then define the interpolated displacements over the element using linear shape functions as ∀𝑥 ∈ [0, 𝐿 𝑒 ]:

⎡ ⎢ ⎢ ⎣ 𝑢 𝑒 (𝑥, 𝑡) 𝑤 𝑒 (𝑥, 𝑡) 𝜃 𝑒 (𝑥, 𝑡) ⎤ ⎥ ⎥ ⎦ = 𝐍(𝑥)𝐪 𝑒 (𝑡), (12)
where the matrix of shape functions 𝐍 is defined as:

𝐍(𝑥) = ⎡ ⎢ ⎢ ⎣ 𝑁 1 (𝑥) 0 0 𝑁 2 (𝑥) 0 0 0 𝑁 1 (𝑥) 0 0 𝑁 2 (𝑥) 0 0 0 𝑁 1 (𝑥) 0 0 𝑁 2 (𝑥) ⎤ ⎥ ⎥ ⎦ , ( 13
)
and the linear shape functions 𝑁 1 and 𝑁 2 given by:

𝑁 1 (𝑥) = 1 - 𝑥 𝐿 𝑒 , 𝑁 2 (𝑥) = 𝑥 𝐿 𝑒 . ( 14
)
Using this interpolation, one can compute the spatial derivatives appearing in the definition of the strain measure in Eq. ( 4):

⎡ ⎢ ⎢ ⎣ 𝑢 𝑒′ (𝑥, 𝑡) 𝑤 𝑒′ (𝑥, 𝑡) 𝜃 𝑒′ (𝑥, 𝑡) ⎤ ⎥ ⎥ ⎦ = 1 𝐿 𝑒 ⎡ ⎢ ⎢ ⎣ 𝑢 2 (𝑡) -𝑢 1 (𝑡) 𝑤 2 (𝑡) -𝑤 1 (𝑡) 𝜃 2 (𝑡) -𝜃 1 (𝑡) ⎤ ⎥ ⎥ ⎦ , ( 15
)
which do not depend on 𝑥 because of the linearity of the shape functions. Substituting the previous relations into Eq. ( 4) leads to the following expressions for the discretized strains over a given element:

𝑒 𝑒 = ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) cos 𝜃 𝑒 + ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) sin 𝜃 𝑒 -1, 𝛾 𝑒 = ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) cos 𝜃 𝑒 - ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) sin 𝜃 𝑒 , ( 16
)
𝜅 𝑒 = 𝜃 2 -𝜃 1 𝐿 𝑒 . Note that at this point the function 𝜃 𝑒 (𝑥, 𝑡) is also discretized and should be considered as:

𝜃 𝑒 (𝑥, 𝑡) = 𝜃 1 (𝑡)𝑁 1 (𝑥) + 𝜃 2 (𝑡)𝑁 2 (𝑥), (17)
which is not substituted into the equations here in order to increase their readability.

The strains in Eq. ( 16) are grouped into the elementary strain vector ε𝑒 defined by the following:

ε𝑒 = [ 𝑒 𝑒 𝛾 𝑒 𝜅 𝑒 ] T . ( 18
)
Then, taking the differential of Eq. ( 16), the variation of the elementary strain vector 𝛿 ε𝑒 can be written:

𝛿 ε𝑒 = 𝐁 𝑒 (𝑥, 𝐪 𝑒 )𝛿𝐪 𝑒 , ( 19
)
with 𝛿𝐪 𝑒 the variation of the nodal displacements and:

𝐁 𝑒 = 1 𝐿 𝑒 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -cos 𝜃 𝑒 sin 𝜃 𝑒 0 -sin 𝜃 𝑒 -cos 𝜃 𝑒 0 𝑁 1 (𝑥)𝐿 𝑒 𝛾 𝑒 -𝑁 1 (𝑥)𝐿 𝑒 (1 + 𝑒 𝑒 ) -1 cos 𝜃 𝑒 -sin 𝜃 𝑒 0 sin 𝜃 𝑒 cos 𝜃 𝑒 0 𝑁 2 (𝑥)𝐿 𝑒 𝛾 𝑒 -𝑁 2 (𝑥)𝐿 𝑒 (1 + 𝑒 𝑒 ) 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ T . ( 20
)
Now that the discretization is defined, it can be substituted into the principle of virtual work of Eq. ( 7) in order to obtain the finite element equations of motion. Beginning with the inertial work component, the expression of Eq. ( 8) over a single element of length 𝐿 𝑒 is first rewritten as:

𝛿𝑊 𝑒 𝑎 = ∫ 𝑉 𝑒 𝜌 ü𝑒 ⋅ 𝛿𝐮 𝑒 d𝑉 = ∫ 𝐿 𝑒 0 𝜌 [ 𝑢 𝑒 (𝑥) 𝑤 𝑒 (𝑥) 𝜃 𝑒 (𝑥) ] 𝐉 ⎡ ⎢ ⎢ ⎣ ü𝑒 (𝑥) ẅ𝑒 (𝑥) θ𝑒 (𝑥) ⎤ ⎥ ⎥ ⎦ d𝑥, (21)
with 𝐉 defined as follows:

𝐉 = ⎡ ⎢ ⎢ ⎣ 𝐴 0 0 0 𝐴 0 0 0 𝐼 ⎤ ⎥ ⎥ ⎦ . ( 22
)
Introducing the discretization leads to the following expression for the virtual work of the inertia forces:

𝛿𝑊 𝑒 𝑎 = ∫ 𝐿 𝑒 0 𝜌(𝐍𝛿𝐪 𝑒 ) T 𝐉𝐍 q𝑒 d𝑥 = (𝛿𝐪 𝑒 ) T ( ∫ 𝐿 𝑒 0 𝜌𝐍 T 𝐉𝐍 d𝑥 ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝐌 𝑒 q𝑒 . ( 23
)
From Eq. ( 23), the elementary mass matrix is derived by evaluating the integral in parentheses, yielding:

𝐌 𝑒 = 𝜌𝐿 𝑒 6 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2𝐴 0 0 𝐴 0 0 0 2𝐴 0 0 𝐴 0 0 0 2𝐼 0 0 𝐼 𝐴 0 0 2𝐴 0 0 0 𝐴 0 0 2𝐴 0 0 0 𝐼 0 0 2𝐼 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 24
)
Following with the virtual work of the internal forces, the expression of Eq. ( 9) is first rewritten as follows:

𝛿𝑊 𝑒 𝑖 = -∫ 𝑉 𝑒 𝐒 ∶ 𝛿𝐄 d𝑉 = -∫ 𝐿 𝑒 0 (𝑁𝛿𝑒 + 𝑇 𝛿𝛾 + 𝑀𝛿𝜅) d𝑥 = -∫ 𝐿 𝑒 0 (𝛿 ε𝑒 ) T Ŝ𝑒 d𝑥. (25)
Introducing the discretization leads to the following expression for the virtual work of the internal forces over a given element:

𝛿𝑊 𝑒 𝑖 = -∫ 𝐿 𝑒 0 (𝐁 𝑒 (𝑥, 𝐪 𝑒 )𝛿𝐪 𝑒 ) T Ŝ𝑒 d𝑥 = -(𝛿𝐪 𝑒 ) T ∫ 𝐿 𝑒 0 (𝐁 𝑒 (𝑥, 𝐪 𝑒 )) T Ŝ𝑒 d𝑥 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝐟 𝑒 int . ( 26
)
The elementary internal force vector is calculated by evaluating the following integral:

𝐟 𝑒 int (𝐪 𝑒 ) = ∫ 𝐿 𝑒 0 (𝐁 𝑒 ) T Ŝ𝑒 d𝑥. (27)
Eq. ( 27) is then readily evaluated using ( 6) and [START_REF] Virgin | Vibration isolation using extreme geometric nonlinearity[END_REF] to yield the elementary internal force vector 𝐟 𝑒 int . Note that the integral is evaluated using a reduced Gauss integration (here with only a single point at 𝑥 = 𝐿 𝑒 2 ) to avoid the phenomenon of shear locking [START_REF] Bathe | Finite Element Procedures[END_REF]. The elementary internal force vector 𝐟 𝑒 int is then given by:

𝐟 𝑒 int = 𝐸𝐴 ē ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -cos θ -sin θ γ 𝐿 𝑒 2 cos θ sin θ γ 𝐿 𝑒 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 𝑘𝐺𝐴γ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ sin θ -cos θ -𝐿 𝑒 2 (1 + ē) -sin θ cos θ -𝐿 𝑒 2 (1 + ē) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 𝐸𝐼 κ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 -1 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 28
)
where • here denotes a quantity evaluated at 𝑥 = 𝐿 𝑒 ∕2, explicitly given by:

θ = 𝜃 1 + 𝜃 2 2 , ( 29
)
ē = 𝑒 𝑒 (𝐿 𝑒 ∕2) = ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) cos θ + ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) sin θ -1, (30)
γ = 𝛾 𝑒 (𝐿 𝑒 ∕2) = ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) cos θ - ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) sin θ, ( 31
)
κ = 𝜅 𝑒 (𝐿 𝑒 ∕2) = 𝜃 2 -𝜃 1 𝐿 𝑒 .
Finally, we consider the virtual work of the external forces in Eq. ( 10) over a given element. First, it is rewritten in the form:

𝛿𝑊 𝑒 𝑒 = ∫ 𝜕𝑉 𝑒 (𝛿𝐮) T 𝐟 d𝐴 + ∫ 𝑉 𝑒 (𝛿𝐮) T 𝐛 d𝑉 = ∫ 𝐿 𝑒 0 ( 𝛿𝑢 𝛿𝑤 𝛿𝜃 ) ⎛ ⎜ ⎜ ⎝ 𝑛 𝑝 𝑞 ⎞ ⎟ ⎟ ⎠ d𝑥, (32)
where 𝑛(𝑥, 𝑡), 𝑝(𝑥, 𝑡) and 𝑞(𝑥, 𝑡) respectively correspond to the distributed axial force, transverse force and bending moment.

Introducing the discretization leads to the following expression for the virtual work of the external forces over a given element:

𝛿𝑊 𝑒 𝑒 = (𝛿𝐪 𝑒 ) T ⎛ ⎜ ⎜ ⎝ ∫ 𝐿 𝑒 0 𝐍 T ⎛ ⎜ ⎜ ⎝ 𝑛 𝑝 𝑞 ⎞ ⎟ ⎟ ⎠ d𝑥 ⎞ ⎟ ⎟ ⎠ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝐟 𝑒 ext , ( 33
)
from which the elementary external force vector 𝐟 𝑒 ext can be derived by evaluating the above integral in parentheses. Gathering the expression for the various virtual works over an element in Eqs. [START_REF] Ben Shaya | Automatic locking of a parametrically resonating, base-excited, nonlinear beam[END_REF], [START_REF] Mahmoodi | Non-linear vibrations and frequency response analysis of piezoelectrically driven microcantilevers[END_REF] and [START_REF] Nayfeh | Linear and Nonlinear Structural Mechanics[END_REF] and summing them over all elements results in the assembly procedure for the elementary matrices and force vectors. If the undeformed structure is curved, one must perform the assembly in the global frame (𝐞 𝑋 , 𝐞 𝑌 ) using a local rotation matrix 𝐓 𝑒 (𝜓) applied to all elementary quantities 𝐪 𝑒 , 𝐌 𝑒 , 𝐟 𝑒 int and 𝐟 𝑒 ext , as explained in Appendix A. After assembly, the discretized nonlinear finite element model falls into the following general form:

𝐌 q + 𝐃 q + 𝐟 int (𝐪) = 𝐟 ext , ( 34
)
where 𝐪 is the column vector of the degrees of freedom expressed in the global frame of the structure, containing all variables 𝑈 𝑥𝑖 , 𝑈 𝑦𝑖 and 𝜃 𝑖 for all the 𝑁 ∈ N nodes, 𝑖 ∈ [1, 𝑁], where the displacement vector of the centerline at position 𝑥 𝑖 is 𝐔(𝑥 𝑖 , 𝑡) = 𝑈 𝑥𝑖 𝐞 𝑋 + 𝑈 𝑦𝑖 𝐞 𝑌 . 𝐌 is the mass matrix of size 3𝑁 × 3𝑁, 𝐟 int (𝐪) is the internal force vector and 𝐟 ext is the vector of externally-applied forces, both of size 3𝑁 × 1. Note that the geometrical nonlinearities are housed within 𝐟 int (𝐪), appearing as sine and cosine functions of the rotation degrees of freedom 𝜃 𝑖 , as shown in Eq. [START_REF] Yang | On the efficiency of piezoelectric energy harvesters[END_REF]. Moreover, in some cases 𝐟 ext can also depend on 𝐪, as in the case of external forces that depend on the deformed geometry of the structure, like follower forces. This case is not explicitly treated in the present study. At this stage, we also define the stiffness matrix of the problem as the tangent stiffness at the equilibrium configuration (𝐪 = 𝟎):

𝐊 = 𝜕𝐟 int 𝜕𝐪 | | | |𝐪=𝟎 . ( 35
)
Note that in Eq. ( 34), a damping matrix 𝐃 has been introduced a posteriori after the FEM discretization. In the test cases that follow, a linear Rayleigh damping model will be considered, where the damping matrix is proportional to the mass matrix:

𝐃 = 𝛼𝐌, ( 36
)
with 𝛼 ∈ R. This kind of damping was chosen here since it is the simplest damping condition to be implemented in the present finite element context and also because this study focuses primarily on the conservative nonlinear modes, as explained in Section 2.3.3.

A full Rayleigh damping model (i.e. 𝐃 = 𝛼𝐊 + 𝛽𝐊 with 𝛽 ∈ R) could also have been selected. However, the stiffness-proportional term creates a modal damping that increases linearly with the frequency, introducing a very heavy damping on the high-frequency modes, which can become overdamped and can subsequently influence the forced responses, likely in a non-physical way (see e.g. [START_REF] Vizzaccaro | Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures[END_REF]). A Kelvin-Voigt damping model could also have been selected, as in [START_REF] Farohki | Geometrically exact extreme vibrations of cantilevers[END_REF][START_REF] Farohki | Experimentally validated geometrically exact model for extreme nonlinear motions of cantilevers[END_REF][START_REF] Amabili | Nonlinear damping in large-amplitude vibrations: modelling and experiments[END_REF]. Since this damping model is coupled to the geometrical nonlinearities, its implementation in our finite element context is not as straightforward as adding a linear damping matrix and is left for upcoming studies. This latter model seems closer to the physics since it perfectly matches measurements in a vacuum [START_REF] Farohki | Experimentally validated geometrically exact model for extreme nonlinear motions of cantilevers[END_REF]. In real, open-air situations, the damping arising from air-structure interactions would be dominant over any structural damping [START_REF] Colin | Very large amplitude vibrations of flexible structures: Experimental identification and validation of a quadratic drag damping model[END_REF], which is even more complex to include in the finite element model and is consequently left for further studies.

Periodic solutions computation and numerical continuation

We turn now to the determination of the solution of Eq. ( 34). Since we are interested in solutions in the frequency domain, we restrict ourselves to periodic solutions, in free oscillation or related to a forced response under harmonic forcing.

Numerically, the harmonic balance method (HBM) is used with the asymptotic numerical method (ANM [START_REF] Damil | A new method to compute perturbed bifurcation: application to the buckling of imperfect elastic structures[END_REF][START_REF] Cochelin | Asymptotic-numerical method and Padé approximations for non-linear elastic structures[END_REF]) to compute the nonlinear modes or the forced response of the nonlinear finite element model. The MANLAB software [START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF][START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], which combines both the HBM and the ANM, is used to carry out the computations of the applications. A key point of the ANM is the quadratic recast [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF][START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF], which is needed to input the equation of motion into MANLAB and which guarantees faster computations and implementations. In this section, each step of the resolution method is presented briefly before moving to the presentation of the test cases.

Quadratic recast for solving with MANLAB

In the remainder of the paper, the ANM is used in conjunction with the HBM to follow periodic solutions (using the MANLAB package). In order to develop an efficient algorithm for the ANM, the HBM is not applied directly onto Eq. [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. Instead, the equation of motion is rewritten in the form of a differential-algebraic system of equations (DAE) with at most quadratic nonlinearities, a process known as the so-called quadratic recast [START_REF] Cochelin | Asymptotic-numerical method and Padé approximations for non-linear elastic structures[END_REF][START_REF] Azrar | An asymptotic-numerical method to compute the post-buckling behavior of elastic plates and shells[END_REF][START_REF] Cochelin | Méthode Asymptotique Numérique, Collections Méthodes Numériques [Eng: The Asymptotic Numerical Method, Collections of Numerical Methods[END_REF]. This step is necessary in order for the resulting DAE system to be of quadratic order, allowing efficient computation of the Fourier coefficients of the nonlinear terms using convolution products. Another advantage of this formalism is that the definition of the Jacobian can be ''hard-coded'' as it only depends on the coefficients of the linear and quadratic parts of the equation. Furthermore, the MANLAB formalism requires that the dynamical system be written in first-order form, in order to be as generic as possible.

The transformation of the FE equation ( 34) into a quadratic DAE involves the definition of new variables, called auxiliary variables. This operation is not unique and several options are possible. Here, we introduce the quadratic recast that leads to the minimum computation time, although it is not the one which introduces the minimum number of auxiliary variables. This is a consequence of the fact that all the auxiliary variables are condensed during computation of the series, such that the number of auxiliary variables has little influence on the computation time. By contrast, the number of product operations of the variables (computing 𝑥 × 𝑦 for variables 𝑥 and 𝑦) increases the computation time and therefore must be minimized, by introducing more auxiliary variables where necessary (see [START_REF] Guillot | A generic and efficient taylor series-based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] for details on the condensation procedure in MANLAB). Outlined below is the quadratic recast that leads to the minimum computation time. Other possible recasts are given in Appendix B along with an analysis of the computation time for each of them.

We first define the following variables:

θ = (𝜃 1 + 𝜃 2 )∕2, (37a)
𝑢 𝑝 = 1 2𝐿 𝑒 (𝑢 2 -𝑢 1 ), ( 37b
)
𝑤 𝑝 = 1 2𝐿 𝑒 (𝑤 2 -𝑤 1 ), ( 37c
)
𝜃 𝑝 = 1 2𝐿 𝑒 (𝜃 2 -𝜃 1 ). (37d)
In addition, the sine and cosine variables 𝑐 and 𝑠 are introduced, defined as:

𝑐 = cos θ, 𝑠 = sin θ. ( 38
)
To render their evaluation quadratic, they are redefined in terms of differentials, given by:

d𝑐 = -𝑠 d θ, d𝑠 = 𝑐 d θ, (39)
which are quadratic w.r.t. the variables 𝑠, 𝑐, d𝑠, d𝑐 and d θ. In order to have Eqs. [START_REF] Reissner | On finite deformations of space-curved beams[END_REF] fully equivalent to Eqs. [START_REF] Reissner | On one-dimensional finite-strain beam theory: The plane problem[END_REF], initial conditions must be added, such as:

𝑐(0) = cos[ θ(0)], 𝑠(0) = sin[ θ(0)], (40)
as explained in [START_REF] Guillot | A generic and efficient taylor series-based continuation method using a quadratic recast of smooth nonlinear systems[END_REF][START_REF] Karkar | A high-order,purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF][START_REF] Guillot | A taylor series-based continuation method for solutions of dynamical systems[END_REF] and detailed on a simple example in Appendix D. Notice that this process is hard-coded in the MANLAB software.

Next, the following additional auxiliary variables are introduced:

𝑒 = ( 1 + 𝑢 𝑝 ) 𝑐 + 𝑤 𝑝 𝑠 -1, (41a)
𝛾 = 𝑤 𝑝 𝑐 - ( 1 + 𝑢 𝑝 ) 𝑠, ( 41b
)
𝐹 𝑥 = 𝐸𝐴𝑒𝑐 -𝑘𝐺𝐴𝛾𝑠, ( 41c
)
𝐹 𝑦 = 𝐸𝐴𝑒𝑠 + 𝑘𝐺𝐴𝛾𝑐, ( 41d
)
𝑀 = 𝐸𝐼𝜃 𝑝 , ( 41e
)
𝑇 2 = 𝐸𝐴𝑒𝛾 -𝑘𝐺𝐴𝛾(1 + 𝑒). ( 41f
)
Note that all of the additional auxiliary variables ( θ, 𝑢 𝑝 , 𝑤 𝑝 , 𝜃 𝑝 , 𝑒, 𝛾, 𝐹 𝑥 , 𝐹 𝑦 , 𝑀 and 𝑇 2 ) are quadratic w.r.t. the variables 𝑢 1 , 𝑢 2 , 𝑤 1 , 𝑤 2 , 𝜃 1 , 𝜃 2 , 𝑐 and 𝑠. Having defined these additional variables, the elementary internal force vector is then rewritten in a quadratic form as a function of the primary and auxiliary variables, as shown below in Eq. [START_REF] Ibrahimbegović | On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements[END_REF].

To summarize, the full quadratic first order DAE is given explicitly for a single finite element. For each element, the following variables are present:

1. Twelve primary variables:

-Six positions:

𝐪 𝑒 = [ 𝑢 1 𝑤 1 𝜃 1 𝑢 2 𝑤 2 𝜃 2 ] T , (42)
-Six velocities:

𝐕 𝑒 = [ 𝑉 𝑢1 𝑉 𝑤1 𝑉 𝜃1 𝑉 𝑢2 𝑉 𝑤2 𝑉 𝜃2 ] T , (43)
2. Twelve auxiliary variables:

𝐗 𝑒 aux = [ θ 𝑢 𝑝 𝑤 𝑝 𝜃 𝑝 𝑐 𝑠 𝑒 𝛾 𝐹 𝑥 𝐹 𝑦 𝑀 𝑇 2 ] T , (44)
yielding the following quadratic first order DAE:

q𝑒 = 𝐕 𝑒 , ( 45a
)
𝐌 𝑒 V𝑒 = 𝐟 𝑒 ext -𝐃 𝑒 𝐕 𝑒 - ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -𝐹 𝑥 -𝐹 𝑦 -𝑀 + 𝑇 2 𝐿 𝑒 ∕2 𝐹 𝑥 𝐹 𝑦 𝑀 + 𝑇 2 𝐿 𝑒 ∕2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟ 𝐟 𝑒 int , ( 45b
) 0 = θ -(𝜃 1 + 𝜃 2 )∕2, (45c)
0 = 𝑢 𝑝 - 1 2𝐿 𝑒 (𝑢 2 -𝑢 1 ), (45d)
0 = 𝑤 𝑝 - 1 2𝐿 𝑒 (𝑤 2 -𝑤 1 ), (45e)
0 = 𝜃 𝑝 - 1 2𝐿 𝑒 (𝜃 2 -𝜃 1 ), (45f)
0 = d𝑐 + 𝑠d θ, 𝑐 = cos θ, (45g) 0 = d𝑠 -𝑐d θ, 𝑠 = sin θ, (45h) 0 = 𝑒 - ( 1 + 𝑢 𝑝 ) 𝑐 -𝑤 𝑝 𝑠 + 1, (45i) 0 = 𝛾 -𝑤 𝑝 𝑐 + ( 1 + 𝑢 𝑝 ) 𝑠, (45j) 0 = 𝐹 𝑥 -(𝐸𝐴𝑒𝑐 -𝑘𝐺𝐴𝛾𝑠), (45k) 0 = 𝐹 𝑦 -(𝐸𝐴𝑒𝑠 + 𝑘𝐺𝐴𝛾𝑐), (45l) 0 = 𝑀 -𝐸𝐼𝜃 𝑝 , (45m) 0 = 𝑇 2 -(𝐸𝐴𝑒𝛾 -𝑘𝐺𝐴𝛾(1 + 𝑒)) . ( 45n
)
In this version of the quadratic recast, there are 12 variables per node and 12 auxiliary variables per element. When considering the full mesh of 𝑁 nodes and 𝑁 e elements, the total number of degrees of freedom of the quadratic first order DAE is then

𝑁 dof = 12𝑁 + 12𝑁 e .

Harmonic balance method

Now that the system of equations is in the form of a quadratic first order DAE, the next step is to apply the harmonic balance method (HBM). Each variable 𝑥 of the 𝑁 dof variables of the DAE (primary and auxiliary variables) is assumed periodic of period 2𝜋∕𝜔, 𝜔 being the angular frequency. The variable 𝑥 is expanded in the form of a truncated Fourier series (up to harmonic 𝐻 ∈ N):

𝑥(𝑡) = 𝑥 0 + 𝐻 ∑ 𝑘=1 (𝑥 𝑐 𝑘 cos 𝑘𝜔𝑡 + 𝑥 𝑠 𝑘 sin 𝑘𝜔𝑡), ( 46
)
where 𝑥 0 , 𝑥 𝑐 𝑘 and 𝑥 𝑠 𝑘 are the Fourier coefficients of the variable 𝑥. In the HBM, the ansatz in Eq. ( 46) is substituted into the quadratic DAE in Eqs. [START_REF] Ibrahimbegović | On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements[END_REF] and the residual is projected orthogonally onto the (truncated) Fourier basis. This leads to an algebraic set of equations that may depend on one or several parameters (excitation frequency, forcing amplitude, etc.). From a general point of view, the system of (quadratic) algebraic equations resulting from the application of the HBM can be written as:

𝐑(𝐗, 𝜔, 𝜆) = 𝐑( X) = 0, (47)
where 𝐗 ∈ R 𝑁 t is the vector containing the Fourier coefficients of all of the variables (primary and auxiliary), of size:

𝑁 t = 𝑁 dof (2𝐻 + 1) = (12𝑁 + 4𝑁 e )(2𝐻 + 1), ( 48
)
𝜆 is a continuation parameter, X = [ 𝐗 T 𝜔 𝜆 ] T ∈ R 𝑁 t +2
and 𝐑 ∶ R 𝑁 t +2 → R 𝑁 t is the function defining the residual of the HBM. The solutions to Eq. ( 47) will be obtained using the asymptotic numerical method (ANM) described hereafter, after describing the two main cases of computations.

Nonlinear modes and forced response

In the applications presented in Section 3, two kinds of solutions will be computed, namely, the periodic response of the system under harmonic forcing (forced response) and the periodic response of the free and undamped system (nonlinear modes).

In the forced response case, the FE equation ( 34) to be solved is as follows:

𝐌 q(𝑡) + 𝐃 q(𝑡) + 𝐟 int [𝐪(𝑡)] = 𝐅 sin 𝛺𝑡, ( 49
)
where 𝐅 is the vector of applied harmonic forcing amplitude and 𝛺 is the angular frequency of the forcing. The computation of the forced response is quite straightforward: the user inputs the DAE Eqs. [START_REF] Ibrahimbegović | On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements[END_REF] and the angular frequency of the periodic solution is chosen equal to that of the forcing. The continuation parameter 𝜆 is also set equal to the value of the frequency, such that 𝜆 = 𝜔 = 𝛺.

The equation 𝜆 = 𝜔 is thus added to the residual Eq. ( 47). The continuation of the periodic solution can then be carried out for several constant forcing amplitudes to reveal the effect of the nonlinearity on the dynamic behavior of the system. In the free response case, the FE equation ( 34) to be solved is initially written as follows:

𝐌 q(𝑡) + 𝐟 int [𝐪(𝑡)] = 𝟎, (50)
where the damping and forcing terms have been canceled. The main differences with the forced response case are that now the system is (a) autonomous and (b) conservative. Due to (a), there is no explicit phase information in the system, so if a periodic orbit of period 2𝜋∕𝜔 is a solution of (50), any periodic orbit identical in amplitude and period, but with a different phase, is likewise a solution. This can lead to numerical instabilities and is traditionally solved with the addition of a phase condition (Eq. (51)b). Furthermore, because of (b), there is an additional variable -the energy of the system -which is implicit and which plays the role of the bifurcation parameter 𝜆. However, because this variable is implicit, the system is overdetermined with the addition of the phase condition. To solve this problem, following [START_REF] Arquier | Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes[END_REF][START_REF] Muñoz-Almaraz | Continuation of periodic orbits in conservative and Hamiltonian systems[END_REF], a fictitious bifurcation parameter 𝜆 is introduced and Eq. ( 50) is replaced with:

{ 𝐌 q(𝑡) + 𝜆 q(𝑡) + 𝐟 int [𝐪(𝑡)] = 0, (a) q𝑖 (0) = 0. (b) (51)
Note that the above system has the form of a classical dissipative autonomous system, well posed for numerical continuation. However, since periodic solutions are computed, the only solution to problem (51) is undamped, and thus the solution computed by the algorithm is 𝜆 ≃ 0 [START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF]. The resulting branches of solution can be plotted in an amplitude/frequency plot in order to visualize the backbone curve of the system.

Asymptotic numerical method and MANLAB

Due to the implicit functions theorem, the solution of Eq. ( 47) with the additional equation (𝜆 = 𝜔 in the forced case or the phase equation ((51)b) in the free conservative case) defines a one-dimensional manifold in the R 𝑁 t +2 space mapped by (𝐗, 𝜆, 𝜔) that can be followed through continuation given an initial point. In this paper, the asymptotic numerical method (ANM), first presented in [START_REF] Damil | A new method to compute perturbed bifurcation: application to the buckling of imperfect elastic structures[END_REF] and derived at length in [START_REF] Cochelin | Asymptotic-numerical method and Padé approximations for non-linear elastic structures[END_REF][START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF][START_REF] Guillot | A taylor series-based continuation method for solutions of dynamical systems[END_REF], is used to carry out the continuation procedure.

The ANM basically relies on a high order predictor. The solution is sought in a parametrized form along with a Taylor series expansion, up to degree 𝑃 :

X(𝑎) = X0 + 𝑎 X1 + 𝑎 2 X2 + 𝑎 3 X3 + ⋯ + 𝑎 𝑃 X𝑃 , ( 52
)
where 𝑎 is a pseudo arc-length and 𝑃 typically chosen equal to 20. In order to obtain the coefficients of this expansion, the ansatz in Eqs. ( 52) is substituted into the HBM residual Eq. ( 47) and terms are gathered relative to the power of the parameter (zero-th order 𝑎 0 = 1, first order 𝑎 1 = 𝑎, second order 𝑎 2 , . . . ). Balancing the coefficients of 𝑎 𝑝 for each order 𝑝 ∈ [0, 𝑃 ] and injecting the definition of 𝑎 results in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]:

1. A nonlinear algebraic system 𝐑( X0 ) = 𝟎 at zero-th order, 2. A cascade of linear systems for the higher orders, which enables the computation of X𝑝 at order 𝑝 ≥ 1 as a function of all of the previous orders X𝑖 , 𝑖 < 𝑝.

Using this procedure, all X𝑝 , 𝑝 ∈ {0, … , 𝑃 } are computed and a branch of solution is obtained as a continuous function of 𝑎. This branch has a radius of convergence 𝑎 m , automatically computed, such that ∀𝑎 < 𝑎 m ‖𝐑( X𝑝 (𝑎), 𝜆 𝑝 (𝑎))‖ < 𝜀, with 𝜀 a tolerance chosen by the user. A numerical continuation process consists then in computing several continuous branches, such that the first point of a given branch X(𝑎 = 0) = X0 is taken as the final point of the previous one, thus identically verifying the zero-th order 𝐑( X0 ) = 𝟎 up to the chosen tolerance 𝜀. The value of X0 for the first branch is computed by a Newton-Raphson algorithm, initiated by an initial guess defined by the user. This latter operation will be called the initialization step in the following section.

In practice, the MANLAB package [START_REF] Guillot | Manlab 4.0: An Interactive Path-Following and Bifurcation Analysis Software[END_REF], implemented in the MATLAB environment, is used to carry out the numerical continuation. In addition, it also offers tools to detect bifurcation points [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF], to compute bifurcated branches and to evaluate the stability of periodic solutions [START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF].

Initialization

Proper initialization of the calculation (i.e., the computation of X0 for the first branch of the continuation) is crucial, as the accuracy of the computation of X0 guarantees the accuracy of all of the branches computed thereafter. Since resonant systems are under investigation here, a preliminary computation of the low frequency eigenmodes is necessary. One solves:

(𝐊 -𝜔 2 𝐌)𝜱 = 𝟎, (53)
to obtain the natural frequencies 𝜔 𝑖 and mode shapes 𝜱 𝑖 for 𝑖 = 1, 2, … When performing forced computations, initializing the Newton-Raphson algorithm with a zero solution far from a resonance works in most cases. Restated more precisely, one prescribes (𝐗 = 𝟎, 𝛺 far from 𝜔 𝑖 ) as the initial guess.

When investigating the nonlinear modes, however, the initialization of a given backbone curve computation is more delicate. Two different options are presented. The first method initializes the algorithm on a linear single-mode solution; namely, for the 𝑖th nonlinear mode 𝐪(𝑡) = 𝛽𝜱 𝑖 cos 𝜔 𝑖 𝑡, with 𝛽 ∈ R small so that the amplitude ‖𝛽𝜱 𝑖 ‖ of the vibration pattern is small enough not to activate the geometrical nonlinearities. The second option, which is the one used in the simulations of the test cases shown in Section 3, is to compute first a harmonically-forced solution branch and then pick a solution (𝐗, 𝜔) close to the 𝑖th phase resonance point, where the displacement is in phase quadrature with the forcing (as done in experimental continuation [START_REF] Denis | Identification of nonlinear modes using phase-locked-loop experimental continuation and normal forms[END_REF]).

Scaling of the equations

To improve the numerical conditioning of the residual 𝐑, it can be useful to scale the parameters such that the equations are rendered dimensionless. Following [START_REF] Neukirch | Comparison of the von Kármán and Kirchhoff models for the post-buckling and vibrations of elastic beams[END_REF], we define the dimensionless variables (denoted by •):

Û𝑥 , Û𝑦 = 𝑈 𝑥 , 𝑈 𝑦 𝐿 , x = 𝑥 𝐿 , t = 𝑡 𝐿 2 √ 𝐸𝐼 𝜌𝐴 , (54a) N, T , q = 𝐿 2 𝐸𝐼 (𝑁, 𝑇 , 𝑞), M = 𝐿 𝐸𝐼 𝑀, n, p = 𝐿 3 𝐸𝐼 (𝑛, 𝑝), ( 54b
)
where 𝐿 is a characteristic length of the structure (e.g. the length of a beam, the diameter of a ring, see examples in Section 3). Rewritten in this way, it can be shown (see Appendix C) that the dimensionless equations of motion of a geometrically exact beam based on Euler-Bernoulli kinematics with uniform cross-section depend on only one parameter, the slenderness ratio: In all computations, 𝑘 = 1 and 𝐺 = 𝐸∕[2(1 + 𝜈)] with 𝜈 = 0.3.

𝜂 = 𝐼∕(𝐴𝐿 2 ). ( 55
)
In the present case of Timoshenko kinematics, the slenderness ratio depends also on the shear parameter 𝜇 defined in Appendix C. The highly flexible structures considered in this article naturally have a very small 𝜂, underlining their slenderness. This indicates that there will be very little effect of shearing of the cross-section even at extreme amplitudes of deformation, meaning that there is virtually no difference between Timoshenko and Euler-Bernoulli kinematics.

Furthermore, introducing dimensionless parameters brings more universality to the results of this article. The dynamical behavior of any structure with the same shape and parameter 𝜂 as the test cases of Section 3 can be quantitatively recovered using the change of variable of Eqs. [START_REF] Mangussi | Resonance in a vibrating beam: a zoo of nonlinear resonance peaks[END_REF].

In practice, for a beam structure with uniform cross-section and homogeneous material, the dimensionless version of system (34) can be obtained simply by setting 𝐸𝐼 = 1, 𝜌𝐴 = 1, 𝜌𝐼 = 𝜂, 𝐿 = 1, 𝐸𝐴 = 1∕𝜂 and 𝑘𝐺𝐴 = 𝑘∕[2(1 + 𝜈)𝜂] with 𝜈 the Poisson's ratio of the material.

Test cases

In this section, the numerical strategy developed in Section 2 is used in the simulation of various flexible beam structures. To highlight the capabilities of the method, different test cases are selected: a simple cantilever beam, a beam clamped at both ends, a circular ring, a truss structure and an idealized model of a tree, as shown in Fig. 3. In each case, the nonlinear modes of the system are computed with Eq. ( 51). The forced responses, however, are only shown for the cantilever beam test case for the sake of efficiency, although they can be easily calculated for each structure.

The nonlinear modes are graphically represented in two ways. Firstly, we show the backbone curve, representing the evolution of the amplitude of the motion as a function of the free oscillation frequency. Precisely, the maximum value of the absolute value over one period of oscillation of the rotation of the cross-section (𝜃(𝑡)) and/or the displacement components of a selected point of the structure ( 𝑈 𝑥 (𝑡), 𝑈 𝑦 (𝑡) ) are depicted. Examples are Fig. 4(a,b,c) for the fundamental nonlinear mode of the cantilever beam. An alternative would have been to show some of the harmonics of these quantities instead of the maximum over a period, in order to investigate the harmonic content of the oscillation. This is addressed with the inclusion of several bar diagrams, mainly in Section 4.

The second graphical representation of the nonlinear modes is their deformed shapes, shown for selected points along the backbone. They are plotted here as several snapshots of the deformed configuration of the structure taken at fixed time over one half period of the periodic motion. Examples are Fig. 4(d,e,f) for the fundamental nonlinear modes of the cantilever beam. The deformed configurations are useful in order to understand more deeply the complex nonlinear dynamics that are observed.

Cantilever beam

We first consider a cantilever beam, with one end fixed and the other free, shown in Fig. 3(a). This configuration allows for very large deformations under small transverse loading conditions, manifesting extreme changes in the shape of the structure. The extreme deformation sometimes even reaches the point of ''bending backwards'', where the free end of the cantilever achieves a displacement behind the opposite fixed end. For this reason, the geometrically exact beam model is a necessity in order to properly capture the geometrically nonlinear behavior. Additionally, the cantilever beam is the simplest example that can be selected in order to validate the numerical strategy proposed in this paper. Many works have been published in the past on nonlinear vibrations of cantilever beams and the interested reader can refer to [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF] for an overview, as well as the recent numerical contributions [START_REF] Shen | Comparison of reduction methods for finite element geometrically nonlinear beam structures[END_REF][START_REF] Vizzaccaro | High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point[END_REF][START_REF] Farokhi | Three-dimensional nonlinear extreme vibrations of cantilevers based on a geometrically exact model[END_REF][START_REF] Farohki | Geometrically exact extreme vibrations of cantilevers[END_REF].

A mesh corresponding to a cantilever beam of length 𝐿 = 1 m, width 𝑏 = 0.05 m and thickness ℎ = 0.001 m, with material properties as shown in Table 1, is generated. Its slenderness ratio is 𝜂 = 8.33 ⋅ 10 -8 , implying a very slender beam. Details on the mesh/system properties, such as the number of discrete elements, are collected in Table 2 for all test cases. For the simulations, as explained in Section 2.3.6, the system's geometry is normalized such that the characteristic length 𝐿 is the length of the beam.

For simulation of the forced response of this system, a uniformly distributed transverse loading 𝑝(𝑥, 𝑡) = 𝑝 0 sin 𝛺𝑡, as defined in Eq. [START_REF] Nayfeh | Linear and Nonlinear Structural Mechanics[END_REF], is prescribed, with 𝑝 0 its amplitude, in [N/m]. It is equivalent to a transverse support motion, in which case the transverse loading becomes 𝑝 0 = -𝜌𝐴 ẅ0 , where ẅ0 is the base acceleration. In the finite element context, this load appears as an external force 𝐟 ext (𝑡) in Eq. ( 49) of the form 𝐟 ext (𝑡) = 𝐅 sin 𝛺𝑡, with the components associated to the transverse displacement in 𝐅 of the form 𝐹 0 [

1 2
1 … 1 1

2

] T and the others being zero. Note that this is valid if the elements are of constant length 𝐿 𝑒 , in which case 𝐹 0 = 𝑝 0 𝐿 𝑒 . In Fig. 4, we indicated its dimensionless value: F0 = 𝐿 2 𝑝 0 𝐿 𝑒 ∕(𝐸𝐼).

The linear damping coefficient 𝛼 of Eq. ( 36) is chosen to be 𝛼 = 2𝜉𝜔 𝑖 with 𝜉 = 0.005 and 𝜔 𝑖 the (linear) natural frequency of the considered mode. Figs. 4-6 depict typical results of our computations for the first three nonlinear modes of the cantilever. Some frequency responses under forced excitation for a given value of the forcing amplitude F0 are shown in blue with the backbone curves added in black. Many interesting nonlinear dynamic phenomena, as were described in Section 1, are uncovered. Beginning with the first nonlinear mode, shown in Fig. 4(a-c), it can be seen that it manifests a slightly hardening behavior, meaning that the resonant frequency 𝛺 increases as the vibration amplitude increases, a characteristic well documented in the literature, see [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF][START_REF] Farokhi | Three-dimensional nonlinear extreme vibrations of cantilevers based on a geometrically exact model[END_REF][START_REF] Pai | Non-linear non-planar oscillations of a cantilever beam under lateral base excitations[END_REF]. Point A represents a point of moderate amplitude of displacement along the backbone curve and point C a point of high amplitude. The solving method introduced in this paper is capable of pushing the cantilever simulation to very extreme amplitudes, as evidenced by the deformed shape at point C. At this point, the cantilever has fully bent backwards, with the free end of the beam vibrating well beyond the location of the fixed end. Note that even though the maximum of the absolute value is shown in Fig. 4(c), the 𝑈 𝑥 displacement is always negative in value as the tip of the beam is always deformed ''backwards'' relative to its initial position.

Point B represents an instance of an internal resonance, wherein some of the energy is transferred to another vibration mode. Typically, in order to characterize the type of internal resonance, one must look at the spatial pattern depicted by each harmonic and compare it to the mode shapes of the structure in order to conclude the type of internal resonance involved, as explained in Section 3.2. In the present case, because this internal resonance is less clear and lacks convergence within the number of retained harmonics, the discussion is left for Section 4.

The second and third nonlinear modes of the cantilever are shown in Figs. 5 and6, respectively. The forced responses and nonlinear modes are shown in 𝜃 for both modes as the rotation of the cross-section 𝜃 gives a good idea of the extent of the beam's deformation. Both the second and third modes of the cantilever beam exhibit a softening trend, wherein 𝛺 decreases with increasing amplitude of vibration, a fact also well documented (see [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF][START_REF] Pai | Non-linear non-planar oscillations of a cantilever beam under lateral base excitations[END_REF]). The amplitude of deformation 𝑈 𝑦 of the beam is less pronounced on the second and third modes than on the first, due to the presence of node points in the deformed shape. Expressed more concretely, due to the particular geometry of the nonlinear mode shapes, a certain maximum rotation amplitude 𝜃 of the cross-section observed at the tip of the beam is obtained for an amplitude 𝑈 𝑦 of transverse displacement, which decreases as the index of the considered mode increases. Since the mechanism of geometrical nonlinearities is linked to the amplitude of the rotations of the cross-sections, the effect of the geometrical nonlinearities is associated with a certain 𝜃 and thus to a 𝑈 𝑦 that decreases with the index of the modes. Another internal resonance point, however, was captured in the simulation of the third mode and appears as an offshoot of the primary backbone curve. This particular internal resonance was not investigated, but a discussion on these internal resonance offshoots appears in Section 5.

Clamped-clamped beam

A second classical type of beam structure is the clamped-clamped (CC) beam. This example is considered here since the mechanism of geometrical nonlinearities is different. For all the other test cases considered in this article, the physical source of geometrical nonlinearities is the large rotation of the cross-sections, but in the case of the CC beam it is the bending/axial coupling (see Fig. 1 and e.g. [START_REF] Lacarbonara | Refined models of elastic beams undergoing large in-plane motions: Theory and experiment[END_REF]). Due to the presence of boundary conditions that restrain axial motion, a bending deformation induces an increase in the axial tension that increases the bending stiffness and creates a strong nonlinear hardening effect. We consider the same beam as in the previous section, but fixed now at its two ends. For the computations, like the cantilever, the CC beam is normalized by its length. The fundamental nonlinear mode is computed with the geometrically exact model. Fig. 7(a) shows the corresponding backbone curve for 𝑈 𝑦 ∕ℎ at the node located at 𝑥 = 0.3𝐿 from the left end of the beam (see Fig. 3(b)). Fig. 7(a) depicts this backbone curve, with a decidedly strong hardening behavior, classical of all flat structures subjected to bending/axial coupling [START_REF] Thomas | Geometrically non-linear flexural vibrations of plates: in-plane boundary conditions and some symmetry properties[END_REF]. This coupling is responsible for the stronger geometrically nonlinear effects than were seen with the cantilever beam. The hardening/softening behavior in particular is significant for much lower amplitudes of transverse displacement when comparing on the order of the thickness ℎ of the CC beam and not of the length of the beam.

The deformed shape of the beam is also shown at two points of the backbone curve. Point F (Fig. 7(b)) represents the classical 1:5 internal resonance (IR) of the first mode with the third mode (see [START_REF] Givois | On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models[END_REF]). The IR is responsible for a strong coupling between these two modes, which is clearly seen in the deformed configuration at this point. The deformed shape exhibits a mix of the first mode shape oscillating at frequency 𝛺 and the third mode shape oscillating at 5𝛺, a pattern visible in Fig. 7(b) as well as in Fig. 13(c). Point G (Fig. 7(c)) is characterized by the deformed shape of the first mode at a high amplitude of displacement, very close to the linear mode shape of mode 1 even for this oscillating frequency (1.5 times higher than the linear natural frequency 𝜔 1 of the structure). Other internal resonances are also computed on the main backbone curve, a 1:9 with mode 5 around 𝛺∕𝜔 1 = 1.56 and a 1:7 with mode 5 around 𝛺∕𝜔 1 = 2.15.

The results of the geometrically exact model and solving method proposed in this paper were compared against the nonlinear von Kármán (VK) analytical model, provided in [START_REF] Givois | On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models[END_REF]. It is solved in [START_REF] Givois | On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models[END_REF] by a modal expansion on the first 8 natural mode shapes of the beam and a numerical continuation of the obtained dynamical system, solved with MANLAB with 𝐻 = 10 harmonics retained in the HBM. The results of this comparison are shown in Fig. 7(a), where the dashed blue line represents the nonlinear VK solution. It can be seen that the geometrically exact model perfectly agrees with the VK model, including the detection of internal resonances during simulation. However, it is noted that the internal resonance ''offshoots'' (known also as internal resonance ''tongues'' [START_REF] Kerschen | Nonlinear normal modes, Part I: A useful framework for the structural dynamicist[END_REF]) are predicted slightly differently in the two models, departing from the central backbone at slightly shifted locations and sometimes in different directions. This can be explained by the different spatial discretizations used in the two models. For the VK model, an analytical modal basis is used with exact natural frequencies, whose ratios to the first natural frequency are: for 𝑖 = 1, … , 5. Considering that the internal resonance stems from the coincidence of one harmonic of the oscillation with the (nonlinear) frequency of a given nonlinear mode (see [START_REF] Givois | On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models[END_REF] for examples), the slight differences between the values of the natural frequencies in the two models explains the discrepancies between the two models. This is illustrated in the inset of Fig. 7(a), in which the shape of the 1:5 IR tongue is computed identically in both models, but with a slight shift in frequency. Fig. 7(a) is also useful to quantify the limitations on the validity of the VK model, since it is based on a truncation of the beam's cross-sectional rotations and is thus valid only for moderate amplitude of vibration. By contrast, the FE geometrically exact model introduced in this paper has no limit of validity. Looking at Fig. 7(a), the backbone computed from the VK model begins to depart from the geometrically exact backbone around 𝑈 𝑦 ≃ 3ℎ.

𝜔 VK 𝑖 ∕𝜔 𝑉 𝐾

Circular ring

Another interesting structure that was tested is the flexible ring. Previous studies have investigated the dynamics of ring structures, notably the studies of Santillan, Virgin and Plaut [START_REF] Virgin | Vibration isolation using extreme geometric nonlinearity[END_REF][START_REF] Santillan | Equilibria and vibration of a heavy pinched loop[END_REF] and, more recently, Lu et al. [START_REF] Lu | Nonlinear vibration isolation via circular ring[END_REF], which were the inspiration for the circular ring structure considered here. The ring of radius 𝑅 = 1 m is discretized into beam elements of circular (radius 𝑟 = 10 mm) cross-section with material properties shown in Table 1 and mesh properties shown in Table 2. Prior to solving, this system is normalized by its characteristic length, which in this case is chosen as the radius 𝑅. The circular ring is clamped on the rightmost node (see Fig. 3(c)).

The resulting dynamical response is shown for the first two nonlinear modes of vibration in Fig. 8. It is readily seen that the first nonlinear mode is a transverse mode, while the second represents an axial mode of vibration. The backbone curves representing the nonlinear modes of the system are traced for the 𝑈 𝑦 ∕𝑅 (for mode 1) and 𝑈 𝑥 ∕𝑅 (for mode 2) displacements at the node directly opposite the fixed node (see Fig. 3(c)). The first nonlinear mode of vibration is obviously softening (Fig. 8(a)). However, an interesting feature is observed in its deformed shape (Fig. 8(b,c)). Typically, if the amplitude of the modal motion were small, related to a linear oscillation, the displacement could be written 𝐮(𝑡) = 𝑎𝜱 cos 𝛺𝑡, where 𝜱 is the eigenvector and 𝑎 ∈ R the amplitude of the motion. Once in the period of oscillation (𝑡 = 0, 2𝜋∕𝛺, 3𝜋∕𝛺, … here), the structure would recover exactly its undeformed shape at this particular time, no matter the extent of the deformation, since 𝐮(0) = 0 ∀𝑎. In other words, when looking at the deformed configurations of the ring in Fig. 8 at points H, I, J and K, the colored lines should pass exactly through the (undeformed) black line. However, at high amplitudes of vibration, the deformed shape of the ring vibrating on its first mode does not return to its initial configuration. Looking at the deformed shape of the ring at point I in Fig. 8(c), the ring does not pass through its initial circular shape, instead maintaining a distorted and slightly elongated oval shape. This is surely linked to a centrifugal effect gaining importance at extreme amplitudes of vibration, due to a ''swinging'' effect of the ring. At lower amplitude (traced at point H), the centrifugal effect is less pronounced, since the ring nearly passes back through its initial circular shape. A second interesting feature is observed in the second mode. Although its backbone curve is globally softening, at low forcing amplitude, the curve is initially of the hardening type, before bending in a softening trend at a higher amplitude. Whereas a switching from softening to hardening behavior as the amplitude increases is classical for slightly curved structures (some examples include shallow shells in [START_REF] Touzé | Non-linear behaviour of free-edge shallow spherical shells: Effect of the geometry[END_REF] and shallow arches in [START_REF] Vizzaccaro | High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point[END_REF][START_REF] Marconi | A higher-order parametric nonlinear reduced-order model for imperfect structures using Neumann expansion[END_REF]), this reversed phenomenon has never been encountered previously, to the knowledge of the authors.

Truss structure

A fourth structure, that of a cantilever-style truss structure, is here considered (see Fig. 3(d)). The analysis of this style of beam was inspired by the work of Guillot, Cochelin and Vergez [START_REF] Guillot | A generic and efficient taylor series-based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] on similar beams. The truss is of length 𝐿 = 1 m and height 𝐿∕10 and is composed of beams of rectangular cross-section of width 𝑏 = 0.1 m and height ℎ = 0.1 m with material properties as shown in Table 1. The truss is divided into four frames of 𝐿∕4, with a diagonal beam connecting the upper left and lower right corners of each frame. The construction of the asymmetrical truss can be seen in Fig. 9, where the solid black lines shown at points L, M and 𝑁 outline the undeformed configuration of the structure. The nodes making up the left edge of the structure are fixed, while all other nodes are free. Like the cantilever, this structure was also normalized by its characteristic length, considered here the length of the structure.

Results of the dynamic simulations are shown for the first two nonlinear modes in Fig. 9. The backbone curves trace the rotation of the cross-section 𝜃 at the lower, right-most node (the bottom right corner of the truss, see Fig. 3(d)). The deformed shapes of the truss are shown at points L and M for the first mode, and 𝑁 for the second mode. The behavior of the truss mirrors the behavior of the cantilever beam, since the first two nonlinear modes of the truss have the same global shape as the first two nonlinear modes of the cantilever. At higher amplitudes, it was observed during simulations (not reproduced here) that some local ''buckling''-like modes, involving single beam members, appeared in the deformed shape. Moreover, one can notice the asymmetrical shape of the deformed shape of the truss, especially visible in Fig. 9(e), for which the maximum (blue) and minimum (green) position of the truss are not symmetrical with respect to the horizontal axis. This is due to the asymmetry in the structure, brought by the four inclined beams that cross the frames.

It is noted that despite the global resemblance between the cantilever beam and the cantilever truss structure, the trend of the backbone curves is reversed: mode 1 of the truss manifests a softening nonlinear behavior, while it is hardening for the cantilever beam; mode 2 of the truss is hardening, whereas it is softening for the cantilever beam. Furthermore, the convergence in harmonics of 𝜃 at point N, shown in Fig. 9(f), is in order to confirm that, even at the relatively high deformation shown at point N, the solution remains converged with only 10 harmonics included Section 4 for more on convergence studies) (see Fig. 10).

Symmetrical tree model

A fifth and final structure is considered in order to confirm that the numerical method outlined in this paper for the nonlinear dynamical simulation of flexible structures is applicable to any finite element structure. For this, a simple finite element model of a tree is introduced. The tree is represented as a collection of several branching beams connected to a central ''trunk''. The tree simulations were inspired by the work of Kovačić et al. [START_REF] Kovacic | Tree-like structures as hierarchical coupled oscillators[END_REF][START_REF] Kovacic | Sympodial tree-like structures: from small to large-amplitude vibrations[END_REF] on the dynamic analysis of tree-like structures based on beam models. The work of Kovačić et al. takes into account in the model many biological inspirations from real trees, such as the diameter, length and mass ratios of branches at each branching point, slenderness coefficient and branching angle. Some of these features are taken into account here to obtain the structure of Fig. 3(e), in particular constant branching angle 𝛼 = 𝜋∕9 and length ratio 𝑙 2 ∕𝑙 1 = 𝑙 3 ∕𝑙 2 = 0.79, as used in [START_REF] Kovacic | Sympodial tree-like structures: from small to large-amplitude vibrations[END_REF]. The ''trunk'' of the tree is of length 𝑙 1 = 0.5 m. The tree structure is constructed of beam elements with constant material properties (see Table 1) and circular cross-section with constant diameter 𝑑 = 0.1 m, simplifying the pattern of decreasing cross-section area of each branch used by Kovačić et al. Unlike the previous four test cases, this model was not normalized by any characteristic length for the simulations. Finally, the bottom-most node is considered fixed in this model, to represent the trunk of the tree being fixed to the ground.

The backbone curves representing the nonlinear modes are traced for the first three nonlinear modes in Fig. 10 at the tip of the rightmost branch, considered as a point that experiences extreme amplitudes of vibration. The deformed shapes of the tree on its first three nonlinear modes are shown at both low and high amplitudes of vibration. Interestingly, the deformed configurations of the first three nonlinear modes of the tree share some similarities with the cantilever beam. The first mode shape of the tree closely resembles the first mode shape of the cantilever. Similarly, the two primary branches of the tree resemble the first mode shape of the cantilever when the tree vibrates on its second mode. The backbone curves of these two modes also manifest hardening behaviors, the same trend as the cantilever's first mode.

By contrast, the third nonlinear mode of the tree mirrors the second nonlinear mode of the cantilever. The deformed configuration of the tree's third mode reveals a set of nodes at the secondary branching points where little to no vibration occurs, which is very similar to the node that appears on the cantilever's second mode. The backbone curve of the tree's third mode also exhibits a softening trend, the same as the cantilever's second mode, furthering the similarities between the two structures.

Convergence studies

An investigation into the convergence of the solutions is here presented based on the harmonic spectra of the variables at various points throughout the simulations. The convergence is presented for the primary variables 𝑈 𝑥 , 𝑈 𝑦 and 𝜃 and the auxiliary variables defined in Eqs. [START_REF] Farokhi | Three-dimensional nonlinear extreme vibrations of cantilevers based on a geometrically exact model[END_REF], [START_REF] Reissner | On one-dimensional finite-strain beam theory: The plane problem[END_REF] and [START_REF] Meier | Geometrically exact finite element formulations for slender beams: Kirchhoff-Love theory versus Simo-Reissner theory[END_REF]. For the solution to be considered accurate, the amplitude of the higher harmonic terms should converge to zero within 𝐻, the number of harmonics retained during the HBM in Eq. ( 46). Table 2 tabulates the number of retained harmonics in each test case.

Fig. 11 depicts the harmonic spectra of Û𝑥 , Û𝑦 and 𝜃 and the auxiliary variables at two points, points A and B, on the backbone of the cantilever beam's first nonlinear mode (see Fig. 4). Due to the geometrical symmetry of the structure in the transverse direction (𝐞 𝑌 ), the harmonic content of the transverse displacement Û𝑦 and the rotation 𝜃 is odd, whereas that of the axial displacement Û𝑥 is even. The spectrum at point A demonstrates an example of a converged solution. In this simulation, 𝐻 = 20 and it is easily seen that the solution has converged within approximately 5 harmonics.

Also in Fig. 11, it is noted that the convergence of the auxiliary variables closely mirrors the convergence of 𝜃, which makes sense as the critical auxiliary variables are those that govern the rotation, 𝑐 and 𝑠, which are based on 𝜃 as shown in Eq. [START_REF] Reissner | On one-dimensional finite-strain beam theory: The plane problem[END_REF]. 𝑈 𝑥 and 𝑈 𝑦 converge ''faster'' than 𝜃, meaning fewer harmonics are needed. This is confirmed in Fig. 12, showing the convergence in 𝐻 at point L on the backbone calculation of the first mode of the truss structure. For this reason, a general and simple rule for checking convergence of a given computation would be to focus on the convergence of 𝜃. The convergence in 𝐻 is shown for the second nonlinear mode of the truss in Fig. 9(f), where even at the higher amplitude of vibration, point 𝑁 on the backbone, the simulation has converged within 𝐻 = 10. Fig. 13 illustrates the harmonic convergence of 𝜃 for two additional points, point C on the backbone of the cantilever's first nonlinear mode, at very high amplitude, and point F, the location of the first internal resonance detected on the backbone of the clamped-clamped beam's first nonlinear mode. Point B's 𝜃 convergence is recopied for comparison.

Throughout all of the calculations, the harmonic content increases as a function of the amplitude of the motion, since the oscillations become more complex and require more harmonics to be properly reconstructed. An example of this trend is visible when comparing the 𝜃 convergence at point A (Fig. 11) on the cantilever's first mode backbone and point C (Fig. 13) further up the backbone: point A converges more rapidly than point C. Beyond point C, convergence is no longer assured within 𝐻 = 20. Care must be taken to ensure that each solution is properly converged. For this reason, the simulations presented in this paper are only presented up to the point where convergence is lost.

A final point to consider is the convergence of the internal resonance tongues, more demanding than the main backbone curve because of the strong coupling with higher harmonics. This point is addressed in the next section and in Appendix E.

Internal resonances

As already explained and illustrated for instance in [START_REF] Touzé | Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques[END_REF][START_REF] Givois | On the frequency response computation of geometrically nonlinear flat structures using reduced-order finite element models[END_REF], an internal resonance (IR) is observed when there is a coincidence between some harmonics of the oscillation frequency of a given mode and the nonlinear frequency of other nonlinear modes. This results in a strong coupling and a radical change of the topology of the backbone curve, giving rise to offshoots, often called IR ''tongues''. We do not consider here the case in which there exists a particular tuning between the (linear) natural frequencies, in which case the backbone curves are completely changed (see e.g [START_REF] Givois | Backbone curves of coupled cubic oscillators in one-to-one internal resonance: bifurcation scenario, measurements and parameter identification[END_REF][START_REF] Gobat | Backbone curves, Neimark-Sacker boundaries and appearance of quasi-periodicity in nonlinear oscillators: Application to 1:2 internal resonance and frequency combs in MEMS[END_REF] for and 1:2 internal resonances).

From our experience and as shown in [START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF][START_REF] Lewandowski | Solutions with bifurcation points for free vibration of beams: an analytical approach[END_REF], there are three topologies for the connections between the IR tongues and the main backbone curve, as illustrated Fig. 14: 1. In type 1, the tongue is connected to the backbone by several bifurcation points that are merged in the plane (amplitude/frequency). The interested reader can refer to [START_REF] Lewandowski | Solutions with bifurcation points for free vibration of beams: an analytical approach[END_REF], which shows nice plots of the unfolding of such bifucation points. This case appears in particular when there are symmetries in the system, such as in the case of the 1:3 IR between modes 1 and 2 of a hinged-hinged beam (see [START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF][START_REF] Lewandowski | Solutions with bifurcation points for free vibration of beams: an analytical approach[END_REF]). 2. In type 2, there are no bifurcations and, during the continuation, the IR tongue appears smoothly by an increase of the harmonics involved in the IR and a sudden bend of the backbone curve. An example of this type is the 1:5 IR of point F in Fig. 7. 3. In type 3, the IR tongue is obtained after symmetry-breaking bifurcations, usually in the case of even IRs such as a 1:2 (see examples in [START_REF] Guillot | A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems[END_REF]).

In practice, computing the IR tongues is often tedious and one usually prefers to avoid it since the main backbone is often the primary concern. In this context, IR types 1 and 3 are interesting in practice, since a basic continuation driver follows the main backbone and the user has the choice to compute the IR tongue or not, using in the former case a bifurcation detector and a branching algorithm (such as the one implemented in MANLAB [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF]). The case of type 2 IR tongues is more tricky since the user must compute a priori the whole loop of the tongue before returning to the main backbone. In such circumstances, the MANLAB solver gives us the capability to ''jump'' over the birth point of the internal resonance on the backbone, from 1 to 2 in Fig. 14. A tangent prediction is used, followed by several Newton-Raphson corrections to converge back to the main branch.

In the test cases presented in Section 3, many type 1 or 3 IRs were detected but not indicated in the figures. Moreover, several type 2 IRs were also detected and, rather than computing the tongues, we chose to ''jump'' over them using the ''jump'' tool of MANLAB. They appear in the figures only with the bottom of their connection to the main backbone curve, so as not to overload the figures and to focus on the main backbone curve.

An unusual case was discovered for the first nonlinear mode of the cantilever beam around 𝛺∕𝜔 1 = 1.031. At this point a strong internal resonance of type 2 was detected, with a particularly large gap between the two connections of the tongue with the main backbone. Moreover, we noticed the smooth apparition of even harmonics in the spectrum, without symmetry-breaking bifurcation (compare points C and B in Fig. 13). Then, following this IR tongue, it was impossible to obtain proper convergence of the number of harmonics within 𝐻 = 20, as seen for Point B in Fig. 13. Looking at the spatial pattern of the harmonics, it seems that a 1:6:15 IR with modes 2 and 3 would be at hand, with some energy spread over several harmonics (3 to 9, including even, and [START_REF] Su | Strain-based geometrically nonlinear beam formulation for modeling very flexible aircraft[END_REF][START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF][START_REF] Gerstmayr | A 3D finite element method for flexible multibody systems[END_REF]. The cantilever's IR is much less clear than the 1:5 IR, Point F of Fig. 13, for which the pattern of the third mode is clearly visible at the fifth harmonic, the energy being largely concentrated in this harmonic. Complete understanding of this particular internal resonance, not measured in the forced regime in [START_REF] Farokhi | Three-dimensional nonlinear extreme vibrations of cantilevers based on a geometrically exact model[END_REF] likely due to damping of the higher modes, is left for further studies.

Conclusion

This paper presented a novel method for the simulation of highly flexible beam structures, based on geometrically exact beam theory and solved using a numerical strategy rooted entirely in the frequency domain. An analytical derivation of the geometrically exact finite element beam model is outlined in the first parts of Section 2. Section 2.3 introduces the solving method. The solver combines principles of the harmonic balance method (HBM) with a continuation method called the asymptotic numerical method (ANM), based in the frequency domain and automated in the MANLAB package.

In order to demonstrate the capabilities of the proposed method, the forced responses and/or nonlinear modes of several test cases were shown in Section 3, including a discussion of the nonlinear phenomena observed in each system and a convergence analysis with respect to the number of harmonics in Section 4. The model is able to capture various nonlinear phenomena, including the classical frequency-dependence on the amplitude of vibration (i.e. hardening and softening trends), bifurcations and instabilities, energy transfers between different modes of vibration, etc. The results of the test cases highlight the power of the method proposed in this paper for the simulation of flexible structures across a variety of applications.

Many advances in the proposed model are planned for the near future in order to better capture the behavior of the flexible structures at extreme amplitudes of vibration. In particular, two main topics will be addressed. First, a stability analysis will be carried out for the free and forced solution, and, second, the model (currently restricted to two dimensions) will be generalized to three dimensions in order to allow for out-of-plane and torsional motion. Finally, future work will aim to validate the proposed model using experimental results and develop a novel nonlinear reduced-order model to reduce computation time.

The authors would like to note also that the finite element model presented in this work is open access and freely available to download in a user-friendly format along with extensive supporting documentation from GitHub [START_REF] Debeurre | 2D geometrically exact (nonlinear) finite element beam model solved in the frequency domain[END_REF].

CRediT authorship contribution statement

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

) 𝑐 + ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) 𝑠 -1, (B.3) 𝛾 = ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) 𝑐 - ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) 𝑠.
To summarize, the full quadratic DAE for a single finite element is stated explicitly. For each element, there are:

1. Twelve primary variables:

-Six positions:

𝐪 𝑒 = [ 𝑢 1 𝑤 1 𝜃 1 𝑢 2 𝑤 2 𝜃 2 ] T , (B .4)
-Six velocities:

𝐕 𝑒 = [ 𝑉 𝑢1 𝑉 𝑤1 𝑉 𝜃1 𝑉 𝑢2 𝑉 𝑤2 𝑉 𝜃2 ] T , (B.5)
2. Four auxiliary variables:

𝐗 4_aux = [ 𝑐 𝑠 𝑒 𝛾 ] T , (B.6)
so that there are 12 + 4 = 18 variables per element. The quadratic DAE is given by (for a single element):

q𝑒 = 𝐕 𝑒 , (B.7a) 𝐌 𝑒 V𝑒 = 𝐟 𝑒 ext -𝐃 𝑒 𝐕 𝑒 - ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝐸𝐴𝑒 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -𝑐 -𝑠 𝛾𝐿 𝑒 ∕2 𝑐 𝑠 𝛾𝐿 𝑒 ∕2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 𝑘𝐺𝐴𝛾 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑠 -𝑐 -(1 + 𝑒)𝐿 𝑒 ∕2 -𝑠 𝑐 -(1 + 𝑒)𝐿 𝑒 ∕2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 𝐸𝐼 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 (𝜃 1 -𝜃 2 )∕2 0 0 (𝜃 2 -𝜃 1 )∕2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟ 𝐟 𝑒 int , (B.7b) 0 = d𝑐 + 𝑠(d𝜃 1 + d𝜃 2 )∕2, 𝑐 = cos[(𝜃 1 + 𝜃 2 )∕2], (B.7c) 0 = d𝑠 -𝑐(d𝜃 1 + d𝜃 2 )∕2, 𝑠 = sin[(𝜃 1 + 𝜃 2 )∕2], (B.7d) 0 = 𝑒 - ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) 𝑐 - ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) 𝑠 + 1, (B.7e) 0 = 𝛾 - ( 𝑤 2 -𝑤 1 𝐿 𝑒 ) 𝑐 + ( 1 + 𝑢 2 -𝑢 1 𝐿 𝑒 ) 𝑠. (B.7f)

B.2. Quadratic recast with 8 auxiliary variables per element

In this recast, we consider the definition of 8 auxiliary variables for each element. With respect to the quadratic recast of Section 2.3.1, the terms 𝐹 𝑥 , 𝐹 𝑦 , 𝑀 and 𝑇 2 are not defined as auxiliary variables. For each element, there are:

1. Twelve primary variables:

-Six positions:

𝐪 𝑒 = [ 𝑢 1 𝑤 1 𝜃 1 𝑢 2 𝑤 2 𝜃 2 ] T , (B .8)
-Six velocities:

𝐕 𝑒 = [ 𝑉 𝑢1 𝑉 𝑤1 𝑉 𝜃1 𝑉 𝑢2 𝑉 𝑤2 𝑉 𝜃2 ] T , (B.
𝐌 𝑒 V𝑒 = 𝐟 𝑒 ext -𝐃 𝑒 𝐕 𝑒 - ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐸𝐴𝑒 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -𝑐 -𝑠 𝛾𝐿 𝑒 ∕2 𝑐 𝑠 𝛾𝐿 𝑒 ∕2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 𝑘𝐺𝐴𝛾 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑠 -𝑐 -(1 + 𝑒)𝐿 𝑒 ∕2 -𝑠 𝑐 -(1 + 𝑒)𝐿 𝑒 ∕2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + 𝐸𝐼 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 -𝜃 𝑝 0 0 𝜃 𝑝 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (B.11b) 0 = θ -(𝜃 1 + 𝜃 2 )∕2, (B.11c) 0 = 𝑢 𝑝 - 1 2𝐿 𝑒 (𝑢 2 -𝑢 1 ), (B.11d) 0 = 𝑤 𝑝 - 1 2𝐿 𝑒 (𝑤 2 -𝑤 1 ), (B.11e) 0 = 𝜃 𝑝 - 1 2𝐿 𝑒 (𝜃 2 -𝜃 1 ), (B.11f) 0 = d𝑐 + 𝑠d θ, 𝑐 = cos θ, (B.11g) 0 = d𝑠 -𝑐d θ, 𝑠 = sin θ, (B.11h) 0 = 𝑒 - ( 1 + 𝑢 𝑝 ) 𝑐 -𝑤 𝑝 𝑠 + 1, (B.11i) 0 = 𝛾 -𝑤 𝑝 𝑐 + ( 1 + 𝑢 𝑝 ) 𝑠. (B.11j)

B.3. Comparison of computation time

In this paragraph, we illustrate the difference in computation time when using the three different formulations of the quadratic recast (resp. with 4, 8 and 12 auxiliary variables). The computations were carried out for the cantilever beam of Section 3.1, discretized with 20 elements and using 20 harmonics in the HBM development. In each case, 50 ANM branches are computed with an ANM threshold set to 10 -6 . Note that for each formulation, the results are the same and only the computation time differs. Table 3 shows the mean computation time per ANM branch for each of the three quadratic formulations presented in this paper.

It can be seen that the smallest computation time is obtained when using the quadratic formulation using 12 auxiliary variables.

As explained in the main text, what increases the computation time is the number of products of variables and not the addition of auxiliary variables. Looking at the proposed DAE, one finds, for each finite element, 12 products of the variables in Eqs. [START_REF] Ibrahimbegović | On finite element implementation of geometrically nonlinear Reissner's beam theory: three-dimensional curved beam elements[END_REF], 24 products in Eqs. (B.7) and 18 products in Eqs. (B.11). Table 3 confirms that the computation time increases with the number of products.

Appendix C. Dimensionless form of the equations

Following [START_REF] Thomas | Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams[END_REF], the strong form of the beam equations equivalent to the weak form of Eq. ( 7) for a geometrically exact beam model with a Timoshenko kinematics is written: where 𝜂 is defined by Eq. ( 55) and:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜌𝐴 ü = (𝑁
𝜇 = 𝐸𝐼 𝑘𝐺𝐴𝐿 2 = 2(1 + 𝜈)𝜂 𝑘 , (C.3)
is the shear parameter [START_REF] Géradin | Mechanical Vibrations, Theory and Application to Structural Dynamics[END_REF], 𝜈 is the material Poisson's ratio and where the second equation is valid for a homogeneous material with 𝐺 = 𝐸∕[2(1 + 𝜈)]. In the case of Euler-Bernoulli kinematics, 𝛾 = 0, so that 𝜇 does not appear in the equations.

Appendix D. Quadratic recast of transcendental function

Let us consider two time functions 𝑥(𝑡) and 𝑦(𝑡) related by a transcendental function. To better understand the ideas, consider here an exponential function: where the first equation is quadratic with respect to 𝑥, 𝑦, d𝑥 and d𝑦.

To illustrate the HBM treatment, we consider only two harmonics in cosine for 𝑥(𝑡) and 𝑦(𝑡): 𝑥(𝑡) = 𝑥 (0) + 𝑥 (1) cos 𝜔𝑡 + 𝑥 (2) cos 2𝜔𝑡, (D.3) 𝑦(𝑡) = 𝑦 (0) + 𝑦 (1) cos 𝜔𝑡 + 𝑦 (2) 𝑦 (1) = 𝑥 (1) 𝑦 (0) , (D.6) 𝑦 (2) = 𝑥 (2) 𝑦 (0) + 𝑥 (1) 𝑦 (1) ∕4. (D.7)

Because of the differentiation in Eq. (D.2a), the zero-th harmonics do not yield an equation. This is corrected by considering Eq. (D.2b) with the Fourier series expansion: 𝑦 (0) + 𝑦 (1) + 𝑦 (2) = exp ( 𝑥 (0) + 𝑥 (1) + 𝑥 (2) ) . (D.8)

As a preliminary conclusion, the above equation proves that the Fourier coefficients of 𝑥(𝑡) and 𝑦(𝑡) are solutions of an algebraic system that is almost quadratic, with only one of the equations (Eq. (D.8) here) incorporating the transcendental function. In the ANM process, this function is treated in the same manner, by differentiating with respect to the arclength parameter 𝑎. Namely, Eq. (D.8) is rewritten:

𝑦(𝑎) = exp[𝑥(𝑎)],
(D.9)

and the two functions 𝑥 and 𝑦 of 𝑎 are expanded in Taylor series (see Eq. ( 52)): which is evaluated by computing 𝐑( X0 ) of the zero-th order of the first series with a Newton-Raphson.

𝑥(𝑎) = 𝑥 0 + 𝑎𝑥 1 + 𝑎 2
In conclusion, it has been shown with this simple example that replacing the transcendental equation (D.1) with (D.2) leads to a quadratic ANM-HBM process enabling computation of all of the unknowns X𝑝 of the problem. This process is hard-coded in the MANLAB software.

Appendix E. Effect of the number of harmonics

In this appendix, we present a quick convergence study with relation to the number of retained harmonics in the HBM development. The computations are carried out for the first nonlinear mode of the cantilever beam presented in Section 3.1, discretized with 20 elements and using a quadratic recast with 12 auxiliary variables. Fig. 15 depicts the evolution of the rotation at the free end of the beam as a function of angular frequency for 𝐻 = 20, 𝐻 = 40 and 𝐻 = 60 harmonics.

It can be seen that the backbone curve of the first mode is well-approximated with 𝐻 = 20 harmonics so as long as the amplitude of the rotation is lower than 3 radians. At higher amplitudes, the harmonic content must be increased in order to obtain a converged solution.

It is also observed that the location of the internal resonance depends on the number of harmonics retained in the HBM development. The higher the number of retained harmonics, the higher the number of uncovered internal resonances, which may be disadvantageous during computation if the primary interest is in the principal resonance curve of a mode. In addition, in order to obtain a correct solution around the internal resonance, a large number of harmonics is required in the HBM development. Even with 𝐻 = 60, the solution around the internal resonance occurring close to 𝛺∕𝜔 1 ≈ 1.03 is clearly not converged, as can be observed in Fig. 16.

For information purposes, Table 4 shows the mean computation time per ANM branch as a function of the number of retained harmonics 𝐻.

Fig. 1 .

 1 Fig. 1. The two main mechanisms of geometrical nonlinearities in slender structures.

Fig. 2 .

 2 Fig. 2. Timoshenko beam kinematics.

Fig. 3 .

 3 Fig. 3. Sketches of the test cases. The red bullet shows for each test case the selected node and the associated displacement/rotation components where backbone curves are plotted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. First nonlinear mode of the cantilever beam of Fig. 3(a). Top row: forced responses ( F0 = 0.05 and F0 = 0.08) and backbone curve of (a) the cross-section rotation 𝜃, (b) the scaled transverse displacement 𝑈 𝑦 ∕𝐿, (c) the scaled axial displacement 𝑈 𝑥 ∕𝐿 at the tip of the cantilever. Maximum amplitude over one period of oscillation. Bottom row: nonlinear deformed shape at (a) low amplitude (point A), (b) internal resonance (point B) and (c) high amplitude (point C); 20 snapshots per half period. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .Fig. 6 .

 56 Fig. 5. Second nonlinear mode of the cantilever beam. Above: forced responses (𝐹 0 = 0.16 and 𝐹 0 = 0.3) and backbone curve of 𝜃 rotation of the cross-section at free end of beam. Bottom: deformed configuration of the cantilever at point D; 14 snapshots per half period. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .

 7 Fig. 7. First nonlinear mode of a beam fixed on both ends (clamped-clamped beam). Top: backbone curve of 𝑈 𝑦 ∕ℎ displacement at a distance 𝑥 = 0.3 [m] from left end of beam (in blue: von Kármán solution). Center: deformed configuration at internal resonance point F. Bottom: deformed configuration at point G; 20 snapshots per half period. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. First and second nonlinear modes of a clamped circular ring. Top row: (left) mode 1 backbone curve of 𝑈 𝑦 ∶ 𝑅 displacement at node opposite the clamp, (center) deformed configuration of the ring at point H, (right) at point I; 14 snapshots per half period. Bottom row: (left) mode 2 backbone curve of 𝑈 𝑥 ∕𝑅 displacement at node opposite the clamp, (center) deformed configuration at point J, (right) at point K; 20 snapshots per half period. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .

 9 Fig. 9. First and second nonlinear modes of a truss structure. Top row: (left) mode 1 backbone curve of 𝜃 rotation of the cross-section at top right corner of free end, (center) deformed configuration of the truss at point L; 2 snapshots per half period, (right) at point M; 4 snapshots per half period. Bottom row: (left) mode 2 backbone curve of 𝜃 rotation of the cross-section at top right corner of free end, (center) deformed configuration at point N; 4 snapshots per half period, (right) harmonic convergence of 𝜃 at point N. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .

 10 Fig. 10. First, second and third nonlinear modes of a tree-like structure. Top row: (left) mode 1 backbone curve of 𝑈 𝑥 displacement at top rightmost node, (center) deformed configuration of the tree at point O; 2 snapshots per half period, (right) at point P; 4 snapshots per half period. Center row: (left) mode 2 backbone curve of 𝑈 𝑥 displacement at top rightmost node, (center) deformed configuration at point Q, (right) at point R; 2 snapshots per half period. Bottom row: (left) mode 3 backbone curve of 𝑈 𝑥 displacement at top rightmost node, (center) deformed configuration at point S, (right) at point T; 2 snapshots per half period.

Fig. 11 .

 11 Fig. 11. Convergence in number of retained harmonics 𝐻 for 𝑈 𝑥 , 𝑈 𝑦 , 𝜃 and auxiliary variables at point A (top) and internal resonance point B (bottom) of cantilever simulation vibrating on its first nonlinear mode.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Convergence in number of retained harmonics 𝐻 for 𝑈 𝑥 , 𝑈 𝑦 , 𝜃 and auxiliary variables at point L of truss simulation vibrating on its first nonlinear mode.

Fig. 14 .

 14 Fig. 14. Types of connections between the internal resonance tongues and the main backbone curve. B indicates a ''bifurcation'' and SB indicates ''symmetry-breaking''. Dashed lines represent unstable parts of the branches.

Marielle Debeurre :

 Debeurre Software, Investigation, Visualization, Writing -original draft, Writing -review & editing. Aurélien Grolet: Software, Conceptualization, Writing -review & editing, Methodology, Supervision. Bruno Cochelin: Conceptualization, Methodology. Olivier Thomas: Conceptualization, Writing -review & editing, Methodology, Supervision.

  cos 𝜃 -𝑇 sin 𝜃) ′ + 𝑛, (a) 𝜌𝐴 ẅ = (𝑁 sin 𝜃 + 𝑇 cos 𝜃) ′ + 𝑝, (b) 𝜌𝐼 θ = 𝑇 (1 + 𝑒) -𝑁𝛾 + 𝑀 ′ + 𝑞, (c) 𝑁 = 𝐸𝐴𝑒, 𝑇 = 𝑘𝐺𝐴𝛾, 𝑀 = 𝐸𝐼𝜃 ′ . (d) (C.1)Using the dimensionless variable definitions of Eqs.[START_REF] Mangussi | Resonance in a vibrating beam: a zoo of nonlinear resonance peaks[END_REF], one obtains:N cos 𝜃 -T sin 𝜃) ′ + n, (a) ŵ = ( N sin 𝜃 + T cos 𝜃) ′ + p, (b) 𝜂 θ = T (1 + 𝑒) -N𝛾 + M′ + q, (c) N = 𝑒∕𝜂, T = 𝛾∕𝜇, M = 𝜃 ′ , (d) (C.2)

  𝑦(𝑡) = exp[𝑥(𝑡)]. (D.1) Differentiating the above equation and adding an initial condition yields the following equivalent system: d𝑦 = 𝑦 d𝑥, (D.2a) 𝑦(0) = exp[𝑥(0)], (D.2b)

  cos 2𝜔𝑡. (D.4) Introducing the above equations into Eq. (D.2a), with d• considered as a time differentiation d•∕d𝑡 and balancing the zero-th and first two harmonics in a standard HBM process leads to:

𝑦 1 =𝑦 2 =

 12 𝑥 2 , 𝑦(𝑎) = 𝑦 0 + 𝑎𝑦 1 + 𝑎 2 𝑦 2 . (D.10) Then, considering Eq. (D.2a) with d• replaced by a differentiation with respect to 𝑎 (d•∕d𝑎) and balancing each power of 𝑎, one obtains: 𝑦 0 𝑥 1 , (D.11) 𝑥 2 𝑦 0 + 𝑦 1 𝑥 1 ∕2. (D.12) Again, because of the differentiation in Eq. (D.2a), one equation is lacking. This is corrected by considering Eq. (D.2b): 𝑦 0 = exp(𝑥 0 ), (D.13)

Fig. 15 .

 15 Fig. 15. Amplitude of the rotation at the free end of the cantilever beam for the first nonlinear mode as a function of the angular frequency for several numbers of retained harmonic (𝐻 = 20, 𝐻 = 40 and 𝐻 = 60).

Fig. 16 .

 16 Fig. 16. Harmonic content of the solution at the internal resonance occurring close to 𝛺∕𝜔 1 ≈ 1.03 computed with 𝐻 = 60 harmonics.

Table 1

 1 Material properties of flexible beam structures, all structures.

	Structure	Density 𝜌 [kg/m 3 ]	Young's modulus 𝐸 [GPa]
	Cantilever	7800	210
	CC	7800	210
	Ring	4400	104
	Truss	7800	210
	Tree	600	13

Table 2

 2 Number of discrete elements per structure and corresponding number of FE nodes.

	Structure	Mode	Elements 𝑁 e	Nodes 𝑁	Harmonics 𝐻
	Cantilever	1	20	21	20
		2	25	26	20
		3	30	31	20
	CC	1	50	51	10
	Ring	1	30	30	10
		2	40	40	10
	Truss	1	85	78	20
		2	85	78	10
	Tree	1	54	55	20
		2	54	55	20
		3	54	55	10

Table 3

 3 Comparison of computation time for the different formulations of the quadratic recast (example of cantilever beam with 20 elements and 20 harmonics).

	No. aux vars per element	4	8	12
	No. of products	24	18	12
	Mean time per ANM branch [s]	2.65	2.63	1.96

Table 4

 4 Comparison of computation time as a function of the harmonic number (example of cantilever beam with 20 elements, quadratic recast with 12 auxiliary variables).

	No. harmonics 𝐻	20	40	60
	No. unknowns in algebraic system	14 022	27 702	41 382
	Mean time per ANM branch [s]	1.96	10.15	15.23

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 860124. The present paper only reflects the authors' view. The European Commission and its Research Executive Agency (REA) are not responsible for any use that may be made of the information it contains.

Data availability

Data will be made available on request.

Appendix A. Elementary quantities in the global frame

This section specifies the change of basis procedure from the local frames to the global frame prior to the finite element assembly process. The angle 𝜓(𝑥) = (ê 𝑋 , 𝐞 𝑥 ) previously introduced defines the orientation of the local frame ( 𝐞 𝑥 , 𝐞 𝑦 ) of a given finite element in the global frame ( 𝐞 𝑋 , 𝐞 𝑋 ) , as shown in Fig. 2. Following discretization of the beam into finite beam elements, the assembly of the elements into the complete structure must be performed in the global frame ( 𝐞 𝑋 , 𝐞 𝑌 ) in order to properly define the deformation of the elements relative to each other. To this end, the displacement degrees of freedom [𝑢 𝑖 , 𝑤 𝑖 , 𝜃 𝑖 ] T , written in the local frame, must be projected in the global frame to obtain [𝑈 𝑥𝑖 , 𝑈 𝑦𝑖 , 𝜃 𝑖 ] T . This is done by using the change of basis matrix:

such that:

where

T gathers the degrees of freedom in the global frame. Inserting the above equation into Eqs. ( 26), [START_REF] Ben Shaya | Automatic locking of a parametrically resonating, base-excited, nonlinear beam[END_REF], and (33) leads to the transformation of all elementary quantities into the global frame, as follows:

prior to assembly, to obtain Eq. [START_REF] Zienkiewicz | The Finite Element Method[END_REF].

Appendix B. Details on the quadratic recast of the equations of motion

The transformation of the FEM Eq. ( 34) into a quadratic DAE involves the definition of new variables, called auxiliary variables. Note that the quadratic recast is not unique. In this section, two other alternative quadratic recasts are given, with less auxiliary variables than the recast presented in Section 2.3.1.

B.1. Minimal quadratic recast with 4 auxiliary variables per element

In this recast, we consider the definition of 4 auxiliary variables for each element, leading to a minimal number of auxiliary variables. We first add the two sine and cosine variables:

The cosine and sine functions are redefined in terms of differential equations, given by: d𝑐 = -𝑠(d𝜃 1 + d𝜃 2 )∕2, d𝑠 = 𝑐(d𝜃 1 + d𝜃 2 )∕2.

(B.2)

We also define the auxiliary variables 𝑒 and 𝛾 related to the axial and shear strains, respectively: