M Claeys

J-J Sinou

J-P Lambelin

B Alcoverro

Multi-harmonic measurements and numerical simulations of nonlinear vibrations of a beam with non-ideal boundary conditions

This study presents a direct comparison of measured and predicted nonlinear vibrations of a clamped-clamped steel beam with non-ideal boundary conditions. A multi-harmonic comparison of simulations with measurements is performed in the vicinity of the primary resonance. First of all, a nonlinear analytical model of the beam is developed taking into account non-ideal boundary conditions. Three simulation methods are implemented to investigate the nonlinear behavior of the clamped-clamped beam. The method of multiple scales is used to compute an analytical expression of the frequency response which enables an easy updating of the model. Then, two numerical methods, the Harmonic Balance Method and a time-integration method with shooting algorithm, are employed and compared one with each other. The Harmonic Balance Method enables to simulate the vibrational stationary response of a nonlinear system projected on several harmonics. This study then proposes a method to compare numerical simulations with measurements of all these harmonics. A signal analysis tool is developed to extract the system harmonics' frequency responses from the temporal signal of a swept sine experiment. An evolutionary updating algorithm (Covariance Matrix Adaptation Evolution Strategy), coupled with highly selective filters is used to identify both fundamental frequency and harmonic amplitudes in the temporal signal, at every moment. This tool enables to extract the harmonic amplitudes of the output signal as well as the input signal. The input of the Harmonic Balance Method can then be either an ideal mono-harmonic signal or a multi-harmonic experimental signal. Finally, the present work focuses on the comparison of experimental and simulated results. From experimental output harmonics and numerical simulations, it is shown that it is possible to distinguish the nonlinearities of the clamped-clamped beam and the effect of the non-ideal input signal.

Introduction

The frequency response of mechanical structures often exhibit nonlinear behaviors such as dependency of eigenfrequencies and dissipation with input amplitude, discontinuities and hysteresis in the frequency response, multi-harmonic response to a mono-harmonic excitation. The sources of these nonlinearities are well known and various simulation methods have been developed to compute a nonlinear frequency response [START_REF] Nayfeh | Nonlinear Oscillations[END_REF], [START_REF] Kerschen | Nonlinear normal modes, Part 1: A useful framework for the structural dynamicist Mechanical Systems and[END_REF]. Yet the first applications of nonlinear simulation to industrial structures are very recent. In particular Renson and Kerschen [START_REF] Renson | Nonlinear modal analysis of the SmallSat spacecraft Topics in Nonlinear Dynamics[END_REF] applied a time-integration method with a shooting algorithm [START_REF] Peeters | Nonlinear normal modes, Part 2: Toward a pratical computation using numerical continuation techniques[END_REF] to compute the nonlinear vibrations of the SmallSat spacecraft, and focused on internal resonances between linear modes. Sinou [START_REF] Sinou | Non-linear dynamics and contacts of an unbalanced flexible rotor supported on ball bearings Mechanism and Machine[END_REF] applied a Harmonic Balance Method to a nonlinear model of an industrial rotor and computed a multi-harmonic frequency response.

Beams are among the few systems whose nonlinear vibrations have been studied both experimentally and theoretically. Nayfeh and Mook [START_REF] Nayfeh | Nonlinear Oscillations[END_REF] studied nonlinear vibrations of beams under several boundary conditions and showed how these dynamical problems can be treated by the method of multiple scales. Tabaddor [START_REF] Tabaddor | Influence of nonlinear boundary conditions on the single-mode response of a cantilever beam[END_REF] pointed out the importance of non-ideal boundary conditions in the modeling of beams. Malatkar and Nayfeh [START_REF] Malatkar | A parametric identification technique for single-degree-of-freedom weakly nonlinear systems with cubic nonlinearities[END_REF] explained how the method of multiple scales enables an easy identification of the nonlinear parameters. This identification method cannot be used on an industrial complex system, but fits well the needs of this study. Kerschen and al. [8] have written an extensive review of nonlinear system identification methods. The two nonlinear methods previously introduced with their industrial applications have also been used for the study of clamped-clamped beam nonlinear vibrations [START_REF] Ribeiro | Non-linear forced vibrations of thin/thick beams and plates by the finite element and shooting methods Computers and Structures[END_REF] [10] [START_REF] Azrar | A semi-analytical approach to the non-linear dynamic response problem of S-S and C-C beams at large vibration amplitudes. Part II: Multimode approach to the steady state forced periodic response[END_REF].

The main objective of this study is to propose a global strategy based on experiments and simulations for identifying the non-ideal boundary conditions and for predicting the nonlinear vibrations of the system that is subjected to multi-harmonic excitations. This global strategy consists in three steps. First of all, the method of multiple scales is used to update the non-ideal boundary conditions of the nonlinear beam with a mono-harmonic excitation. Secondly, the Harmonic Balance Method is developed to compute easily the multi-harmonic frequency response of the system to multiharmonic excitations. Finally, the third step consists in comparing the multi-harmonic response with experiments. Such a comparison needs an efficient signal processing tool to extract the experimental multi-harmonic frequency response from a temporal signal. This paper proposes and presents a process based on high-selective filtering and fitting to synthetic functions using the Covariance Matrix Adaptation Evolutionary Strategy (CMAES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF].

In addition, the Harmonic Balance Method is compared to the time-integration method with a shooting algorithm [START_REF] Peeters | Nonlinear normal modes, Part 2: Toward a pratical computation using numerical continuation techniques[END_REF], to demonstrate the efficiency of both methods, and complementary experimental measures are performed to substantiate the model developed.

The next section of this article presents the clamped-clamped beam experiment with the method used to compute the experimental multi-harmonic frequency response. Then the beam model with non-ideal boundary conditions is detailed. The fourth section of this article focuses on the updating of the non-ideal boundary conditions by the method of multiple scales. The two numerical methods are explained in the fifth section. Finally, the multi-harmonic comparison between experimental measures and simulations is presented and interpreted. The simulations are made both with an ideal monoharmonic input and with the experimental multi-harmonic input.

Experiments

Experimental setup

The experimental setup is presented in Figure 1. The system studied is a steel beam of dimension 470 x 20 x 5 mm. At both ends, the beam is bonded to a heavy steel block, of dimension 100 x 100 x 85 mm. The whole piece (beam and blocks) is manufactured from a unique bulk piece of steel. Hollows between the beam and the blocks are designed to avoid stress concentration. The blocks are screwed on a large circular aluminum plate, itself screwed on the shaker. The system is instrumented with a strain gage and 3 three-dimensional accelerometers as depicted in Figure 1. Each signal is denoted by the sensor name followed by the direction of the measurement. The entrance signal 1 2 (P1Z+P2Z) is the shaker's feedback control signal. The reference frame of the study is the one of the blocks, where P1 and P2 are fixed. The zero-point of this reference frame coincides with A1 when the system is at rest. The response signal is the acceleration at the middle of the beam, in the reference frame of the blocks, that is A1Z -1 2 (P1Z+P2Z). The horizontal displacement of the blocks is also measured with the signal 1 2 (P2X-P1X). The strain gage is used to evaluate the static constraints in the beam (see Section 2.2). Increasing and decreasing swept sine experiments in the vicinity of the primary resonance is realized.

Experimental results and discussion

Experimental results for three levels of excitation are plotted in Figure 2. The arrows indicate the directions of sweep, which highlight the jump and hysteresis phenomena.

Beyond any nonlinear phenomenon, some variations of the fundamental frequency of the beam (identified at low excitation level) are observed from one experiment to another. A strain gage is placed on the beam (see figure 1) to understand the origin of this phenomenon. The frequency shifts are well correlated with the static value of the strain gage. The static strain evolves a lot during the screwing of the blocks, in a non-reproducible way. It is assumed that the frequency variations come from the position of the blocks that evolve during screwing and from one experiment to another. This movement induces static constraints in the beam. The fundamental frequency of the beam shall be updated for each experiment. The updating of the eigenfrequency is made through an effective half-length of the beam l (235 < l < 275), that takes into account the influence of both the hollows at the ends of the beam and the constraints.

On the contrary, the experimental response shape, which characterizes the nonlinear phenomena appears to be very reproducible from one experiment to another. This visual diagnosis is confirmed by updating the model (see Section 4). The effective half-length l is the only parameter that does not remain stable from one experiment to another.

Multi-harmonic response measurements

The vibrator is piloted with a logarithmic swept sine signal. The frequency evolves slowly enough (0.1 octave/min in the vicinity of the resonance) to consider the response to be stationary. In the vicinity of each time point, both input and output signals are periodic and can thus be developed in Fourier series. The signal analysis tool developed in this study aims at identifying at each time the fundamental frequency and the harmonics' amplitudes of both input and output signals. It is Frequency (Hz) Amplitude of the harmonic response (m) Figure 2: Harmonic response of the beam to 3 levels of excitation signal: in blue 2 m.s -2 , in green 4 m.s -2 and in red 6 m.s -2 . separated in four main steps; the windowing of the temporal signals, the identification of the mean fundamental frequency of each time slot, the band-pass filtering of the signals around the harmonic frequencies and finally the identification of the harmonic amplitudes. The overall process is summed up in figure 3. The windowing step drives the frequency increment of the response curve. The temporal distance between 2 time slots can be directly linked to a frequency increment through the logarithmic speed of sweep. In this work, a constant distance of 0.5s between 2 time slots is chosen. Then the time slot duration is chosen. There should be enough points to perform the identification steps but the corresponding frequency length of the slot should be low enough to consider the harmonic amplitudes to be constant along the sample. In this work, a constant length corresponding to 5 periods of the lowest excitation frequency of the swept sine experiment is chosen. The frequency identification step is detailed in Figure 4. The mean fundamental frequency of each slot is first estimated from the pilot's sweep profile. This prediction is then updated. In order to cancel the higher order harmonic signals, the input signal is filtered by a band-pass filter around the estimated fundamental frequency F est . A band-pass Finite Impulse Response (FIR) filter between 0.9 F est and 1.1 F est is chosen because its effect on the signal phase can be corrected by a simple time translation. The length of the FIR is a compromise between rejection of the higher order harmonics and attenuation in the selected frequency band, a length of 1000 is chosen. In order to get rid of the transient signal introduced by the filter, a larger sample of the input signal is filtered and the original slot is then selected in this filtered sample. Once filtered, the input signal is mono-harmonic; the fundamental frequency can then be identified by fitting a synthetic signal. The synthetic signal chosen is a 4 parameter logarithmic swept sine function ("chirp" function). The four parameters are the constant amplitude A, the frequency at the beginning and the end of the sample F beg , F end and the initial phase Φ 0 .

Band-pass filtering

ϕ chirp (t) = A sin   2π   F end F beg t t end F beg   t + Φ 0   (1)
This identification is treated as a nonlinear optimization problem. The optimization function is the root mean square difference between the sample and the synthetic signal. The four parameters are evaluated, bounded and then optimized using an evolutionary algorithm: the "Covariance Matrix Adaptation -Evolution Strategy" (CMAES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] which is chosen because of its quick convergence and its weak noise sensitivity. The amplitude A is estimated by the maximal absolute value in the sample, the estimated value of both F beg and F end is F est . As the phase Φ 0 is a very sensitive parameter, it is previously estimated by a zero-crossing scheme. The beginning of the synthetic signal is set at the first point after zero crossing. There, the phase Φ 0 is estimated by ± arccos( y A ), where y is the value of the first point, the sign depends on the direction of the zero crossing. Once estimated, the parameters are bounded for the optimization algorithm. A, F beg and F end are bounded to their estimation ±10%, Φ 0 is bounded to its estimation ±10 degrees. At the end of the optimization process, the fundamental frequency F of the sample is defined as the logarithmic mean between F beg and F end .

For each harmonic n computed, the samples of input and output signals are filtered around n.F using the FIR filter described previously. Each filtered sample is then fitted to a two parameter swept sine function, which is the same chirp function presented previously where the beginning and end frequencies are fixed as n.F beg and n.F end . The amplitudes identified are the harmonic amplitudes H n input and H n output . For the output signal, the acceleration amplitudes are translated into displacement by dividing by (2.π.n.F ) 2 . The experimental harmonic curves obtained by this process are plotted in Figure 5. Note that this identification method also enables to extract harmonics' phase information, though this information is not studied in this work.

Other identification methods could have been used (for example methods based on Discrete Fourier Transform or Auto-Regressive formulations [START_REF] Nield | A review of time-frequency methods for structural vibration analysis[END_REF]). The optimization algorithm measures the error between experimental and fitted data and thus ensures the precision of the whole identification process. On synthetic multi-harmonic signals, this identification method provides errors < 0.1% for the frequencies and < 1% for all harmonics amplitudes.

Modeling

The beam equations developed by Nayfeh [START_REF] Nayfeh | Nonlinear transverse vibrations of beams with properties that vary along the length[END_REF] are used and the beam properties are assumed to be constant along the length. The low thickness-length ratio of the beam and the boundary conditions enable to neglect the inertial and curvature nonlinear terms pointed out by Nayfeh ([14], Equation ( 14)). A linear viscous damping is added in order to take into account dissipation in the dynamical differential equation:

ρA ∂ 2 w ∂t 2 + µ ∂w ∂t + EI ∂ 4 w ∂x 4 = ρA a excit cos(Ω t) + T (t) ∂ 2 w ∂x 2 (2)
The variables are defined as follows: ρ is the mass density, A is the cross-sectional area, µ is the linear viscous damping coefficient, E is Young's modulus, I is the cross-sectional moment of inertia, a excit is the acceleration of both clamped ends of the beam, t is time, x is the abscissa, w is the transverse displacement in the reference frame of the clamped ends, u is the axial displacement and T (t) is the tensile force in the beam, invariant along the length, that satisfies the relation:

T (t) = EA ∂u ∂x + 1 2 ∂w ∂x 2 ∀x ∈ [0, l] (3)
As the structure and its excitation are symmetric, the response w is symmetric and only half of the beam shall be studied. The boundary conditions in x = 0 and x = l are detailed in equations ( 4) and ( 5), and illustrated in Figure 6. These non-ideal boundary conditions are due to the deformable block and hollow at each end of the beam. The springs k rot and k bound , the viscous damping coefficient µ and the effective half-length l are unknown parameters that are updated in Section 4. The values of the other physical parameters are given in Table 1. The accelerometer A1 is not modeled as a concentrated mass in order to simplify the modeling. Its influence on the dynamic behavior of the beam is weak (m A1 = 5.7g << m ef f beam = 304g, where m ef f beam is the effective mass of the beam fundamental mode) and is mainly taken into account through the effective half-length parameter l during the model updating. Considering Equations ( 3) and ( 5), we can express T as a function of w and write the transverse dynamical equation [START_REF] Tabaddor | Influence of nonlinear boundary conditions on the single-mode response of a cantilever beam[END_REF], that defines the dynamical nonlinear problem once associated to the boundary conditions [START_REF] Peeters | Nonlinear normal modes, Part 2: Toward a pratical computation using numerical continuation techniques[END_REF].

w(0) = 0, EI ∂ 3 w ∂x 3 (0) = k rot ∂w ∂x (0), ∂w ∂x (l) = 0, ∂ 3 w ∂x 3 (l) = 0 (4)
T (t) = k bound u(0) and u(l) = 0 (5)
ρA ∂ 2 w ∂t 2 + µ ∂w ∂t + EI ∂ 4 w ∂x 4 = ρA a excit cos(Ω t) + EA 2l 1 + EA l k bound -1 l 0 ∂w ∂x 2 dx ∂ 2 w ∂x 2 (6)
The numerical simulation methods, such as the Harmonic Balance Method used in this work, can only solve a set of discrete differential equations. The continuous equation ( 6) is thus projected on the modal basis of the associated homogeneous equation [START_REF] Malatkar | A parametric identification technique for single-degree-of-freedom weakly nonlinear systems with cubic nonlinearities[END_REF].

ρA ∂ 2 w ∂t 2 + EI ∂ 4 w ∂x 4 = 0 (7)
The modal solution of ( 7) is:

Y i (x) = α sin λ i x l + β cos λ i x l + γ sinh λ i x l + +δ cosh λ i x l (8)
with {α, β, γ, δ} solution of:

     0 1 0 1 1 + EIλ 2 i krotl 2 0 1 - EIλ 2 i krotl 2 0 cos(λ i ) -sin(λ i ) cosh(λ i ) sinh(λ i ) -cos(λ i ) sin(λ i ) cosh(λ i ) sinh(λ i )           α β γ δ      =      0 0 0 0      (9)
and λ i is the i th solution of the equation:

1 - EIλ 2 i k rot l 2 cos(λ i ) sinh(λ i ) + 1 + EIλ 2 i k rot l 2 sin(λ i ) cosh(λ i ) = 0 (10)
Y i can be rewritten as: In the experiment, the transverse displacement is measured at the middle of the beam. The modal projections are normalized at this point in order to get the simple relation w (l) = w i . The set of indexes corresponding to the selected modes is named M . As only the first resonance is studied, the set M = {1, 2} is chosen, the shapes of the selected modes are plotted in Figure 7. The projection gives the set of equations [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF].

Y i (x) = -tan(λ i ) sin λ i x l -cos λ i x l -tanh(λ i ) sinh λ i x l + cosh λ i x l (11
d 2 w i dt 2 + µ ρA dw i dt + Ω 2 i w i = Γ i a excit cos(Ωt) + {j,k,m}∈M 3 Γ ijkl . w j w k w m ∀i ∈ M (12)
where:

w i (t) = Y i (l) l l 0 w(x, t).Y i (x) dx; Ω i = EI ρA λ 2 i l 2 ; Γ i = Y i (l) l l 0 Y i (x) dx Γ ijkm = -E Y i (l) 2ρ l 2 Y j (l)Y k (l)Y m (l) 1 + EA l k bound -1 l 0 dY j dx dY k dx dx l 0 dY m dx dY i dx dx
4 Updating via the method of multiple scales

The method of multiple scales, as described by Nayfeh and Mook [START_REF] Nayfeh | Nonlinear Oscillations[END_REF], enables to solve analytically the nonlinear problem (6) in the vicinity of the primary resonance. In this study, this approach is used to update the non-ideal boundary conditions of the beam. As previously explained, the unknown parameters are the effective half-length of the beam l, the viscous damping coefficient µ, the rotational boundary spring k rot and the longitudinal boundary spring k bound .

The variables are non-dimensionalised according to the following scheme:

t = EI ρA λ 2 l 2 t; Ω = ρA EI l 2 λ 2 Ω; x = x l ; r = I A (radius of gyration); w = w r ; ãexcit ( t) = ρAl 4 EIλ 4 r a excit (t); μ = l 2 λ 2 √ ρAEI µ; kbound = l EA k bound ; krot = l 2 EIλ 2 k rot (13)
An arbitrary non-dimensional small parameter ε is introduced, it represents the order of magnitude of the small displacement w, then t and w are expanded as polynomials in ε:

t = t0 + ε t1 + ε 2 t2 + ... and w = ε w1 + ε 2 w2 + ε 3 w3 + ... (14)
As the nonlinear term of the dynamical equation ( 6) has the order of magnitude of w3 which is ε 3 , the damping and the excitation terms are redefined in order to make their order be also ε 3 :

μ = μ ε 2 and ãexcit = ãexcit ε 3 (15)
The method of multiple scales assumes that the various scales ε i are independent so that the dynamical equation ( 6) can be split to get the system [START_REF] Blevins | Formulas for natural frequency and mode shape Krieger publishing company[END_REF].

1 st order ∂ 2 w1 ∂ t2 0 + ∂ 4 w1 ∂ x4 = 0 2 nd order ∂ 2 w2 ∂ t2 0 + ∂ 4 w2 ∂ x4 = -2 ∂ 2 w1 ∂ t0 ∂ t1 3 rd order ∂ 2 w3 ∂ t2 0 + ∂ 4 w3 ∂ x4 = -2 ∂ 2 w1 ∂ t0 ∂ t2 - ∂ 2 w1 ∂ t2 1 -2 ∂ 2 w2 ∂ t0 ∂ t1 - μ ∂ w1 ∂ t0 + ãexcit cos( Ωt ) + 1 + 1 kbound -1 1 0 ∂ w1 ∂ x 2 dx ∂ 2 w1 ∂ x2 (16)
∀i ∈ {1, 2, 3}, wi (0) = 0;

∂ 3 wi ∂ x3 (0) = krot ∂ wi ∂ x (0); wi (1) = 0; ∂ wi ∂ x (1)
= 0 w1 is a solution of the homogeneous problem. As the problem is studied in the vicinity of the first resonance, it is assumed that w1 is collinear to the first modal shape Ỹ1 (x) = Y 1 (x). w1 (x, t0 , t1 , t2 ) = Ỹ1 (x)A 1 ( t1 , t2 )e i t0 + cc (complex conjugate) [START_REF] Cameron | An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems[END_REF] The independence between the scales w i and the uniformity of such an expansion rely on their orthogonality and on a solvability condition. As the problem is self-adjoint (which can be demonstrated following the procedure outlined by Nayfeh [START_REF] Nayfeh | Perturbation Methods[END_REF]), the solvability condition at each order is the elimination of the secular terms. At the second order the solvability condition is ∂A 1 ∂ t1 = 0. Thus w2 is a solution of the homogeneous problem and is orthogonal to w1 . w2 is null. The detuning parameter σ is introduced to measure the distance Ω -Ω 1 . It is defined by:

Ωt = t0 + σ t2 ( 18
)
This definition is equivalent to its dimensional form:

σ = Ω -Ω 1 ε 2 Ω 1 (19)
The solvability condition at the third order, considering the expressions ( 17) and ( 18), gives the equation:

-i dA 1 d t2 -i μ 2 + a excit 2 e iσ t2 -3A 2 1 Ā1 1 + 1 kbound -1   1 0 d Ỹ1 dx 2 dx   2 = 0 (20)
By introducing the amplitude a 1 (such as A 1 = a 1 2 e iβ ), splitting the real and imaginary parts of Equation ( 20) and considering the stationary solution, the relation (21) between the amplitude of vibration a 1 and the frequential detuning σ is obtained.

σ = - 3 4 1 + 1 kbound -1   1 0 d Ỹ1 dx 2 dx   2 a 2 1 ± ãexcit a 1 1 0 Ỹ1 dx 2 - μ 2 2 (21)
In order to compare the simulation results with experiments, Equation ( 21) is rewritten into a dimensional relation ( 22) between the frequency of excitation F and the amplitude of displacement at the middle of the beam w. The corresponding curve is plotted in Figure 8.

F = F 0   1 + B w 2 ± C a excit Ω 0 2 w 2 -D 2   (22)
where

F 0 = 1 2π EI ρA λ 2 1 l 2 ; Ω 0 = 2πF 0 ; B = -3 Al 2 16IY 1 (l) 2 λ 4 1 l 0 dY 1 dx 2 dx 2 1 + EA l k bound -1 C = Y 1 (l) l l 0 Y 1 (x)dx ; D = l 0 Y 1 (x)dx -1 l 3 λ 2 1 Y 1 (l) √ ρAEI µ (23)
In the model presented in Section 3, four parameters have to be experimentally updated: the linear damping coefficient µ, the effective half-length l and the boundary springs k bound and k rot . To this purpose, the following strategy is applied. Equation ( 22) gives a direct link between the frequency response curve and the parameters. The physical terms F 1 , F 2 , F max and w max , defined by Figure 8, are identified in the experimental frequency response curve. The curve parameters F 0 , B, C and D are deducted from relations (24) to (27).

F 0 = F 1 + F 2 2 (24) B = F max -F 0 F 0 w 2 max (25) D = a excit w max Ω 0 2 (26) C = a excit Ω 0 2 w low 2 -D 2 -1 2 F 2 -F 1 2 (27)
The model parameters µ, l, k rot and k bound are then computed using the relations (23). The dependences are the following : C(k rot ), F 0 (k rot , l), B(k rot , l, k bound ), D(k rot , l, µ). Using the relation C(k rot ), k rot is firstly adjusted by a bisection method. l, k bound and µ are then directly computed Table 2: Physical values of the updated parameters.

F F max F 2 F = F 0 (1 + Bw 2 ) w low w F 1 w max F 0
using the relations (23). Table 2 gives the values of four updated parameters for each experiments. It can be noted that only l varies from one experiment to another. The variations of l can be interpreted as variations of the static constraint in the beam from one experiment to another and then correlated to the value of the strain gage during each experiment. When an additional static longitudinal loading δP is added to the system, the fundamental frequency is modified according to the following equation ( [START_REF] Blevins | Formulas for natural frequency and mode shape Krieger publishing company[END_REF], p144):

δF 0 F 0 = 1 2 δP P buckling ( 28
)
where P buckling = π 2 EI l 2 for a clamped-clamped beam experiment ([16], p149). From Equation (23), this fundamental frequency variation δF 0 can be translated into an effective half-length variation δl: The offset for the experimental strain values is chosen to set the mean value of the first experiment to zero. Resonance phenomena appear at 100s (increasing sweep) and 480s (decreasing sweep).

δF 0 F 0 = -2 δl l (29)
The loading δP can be translated into a longitudinal strain δε:

δP = EA δε (30)
Finally, δl is related to the strain δε through the following equation:

δε = -4 P buckling EA δl l (31)
The updated δl are translated into δε and are compared to the experimental data in Figure 9. The experimental data are the strain gage values, low-pass filtered and averaged. Note that in the vicinity of the resonance phenomena, the static value of the gage increases, this is due to the longitudinal movement u which is proportional to w 2 cos 2 (Ωt) and thus have a static component. Away from these resonances, as w 2 is very low, the static value of the gage is not influenced by the dynamical response.

The strain evolves a lot during the experiment, so that it is very difficult to define precisely a constant longitudinal loading. Anyway the updated values of such a loading, obtained via the method of multiple scales, are compatible with the measures. This comparison substantiates the explanations of Section 2.2 on the variability of the fundamental frequency and allows to update l at each experiment.

Numerical Methods

As the method of multiple scales does not allow to simulate a frequency response on a large frequency range, or the response of an industrial finite element model, two numerical methods are employed in this study.

Harmonic Balance Method

The Harmonic Balance Method [START_REF] Sinou | Non-linear dynamics and contacts of an unbalanced flexible rotor supported on ball bearings Mechanism and Machine[END_REF], enables to solve numerically the discrete problem [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], which can be rewritten in the matrix form:

Ẅ + D Ẇ + KW = F excit + F nl (W) ( 32
)
where

W = {w i } i∈M ; D = µ ρA ; K = Diag(Ω 2 i ) i∈M ; F excit = {Γ i a excit cos(Ωt)} i∈M F nl (W) = {j,k,m}∈M 3 Γ ijkm . w j w k w m i∈M (33)
As the excitation term is periodic, it is assumed that the nonlinear dynamical response and the force vector may be approximated by finite Fourier series with Ω as fundamental frequency.

W(t) = B 0 + p k=1 (B k cos(kΩt) + A k sin(kΩt)) (34)
F nl (t) = C 0 + p k=1 (C k cos(kΩt) + S k sin(kΩt)) (35) F excit (t) = p k=1 (C k, excit cos(kΩt) + S k, excit sin(kΩt)) (36) with C 1, excit = {Γ i a excit } i∈M , C k, excit = 0 ∀k = 1 and S k, excit = 0 ∀k ∈ [1, p]
p, the order of the Fourier series is selected on the basis of the number of significant harmonics expected in the nonlinear dynamical response. For our application, p = 5 is chosen. Equations [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] are rewritten in the Fourier basis:

KB 0 = C 0 and K -(kΩ) 2 I d -kΩD -kΩD K -(kΩ) 2 I d A k B k = S k C k + S k,excit C k,excit ∀k ∈ [1, p]
(37) where I d is the identity matrix.

The coefficients C 0 , S k and C k are calculated by applying an Alternate Frequency Time domain method (AFT-method) that was presented by Cameron and Griffin [START_REF] Cameron | An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems[END_REF]. This method is summed up by the following scheme, where DFT means Discrete Fourier Transform:

X = [B 0 A 1 B 1 . . . A p B p ] DF T -1 → W(t) eq(33) → F(t) DF T → [C 0 S 1 C 1 . . . S p C p ] ( 38
)
This link between [C 0 S 1 C 1 . . . S p C p ] and X enables a nonlinear solving of problem (37) which can be written as H(X, Ω) = 0. In order to compute a continuous response curve in the frequency-amplitude domain, a pseudoarclength continuation method is used. The frequency Ω is considered as a variable as well as the Fourier coefficients X. The arc length between two consecutive points (X (i) , Ω (i) ) and (X (i+1) , Ω (i+1) ) is given by:

∆s i+1 = X (i+1) -X (i) 2 + α Ω (i+1) -Ω (i) 2 (39)
α is a numerical coefficient introduced in order to take into account the difference of magnitude between Ω (∼ 10 3 s -1 ) and X (∼ 10 -3 m). This coefficient plays a key role in the numerical convergence of the algorithm, in particular in the vicinity of a singular point. The continuation method is illustrated in Figure 10. It starts with a prediction step. Considering that 3 points are already computed (Y i-2 , Y i-1 and Y i , where Y i = [X (i) , Ω (i) ]), the next point is extrapolated. A second degree Lagrangian polynomial prediction is chosen, so that the prediction point at the abscissa s i+1 is:

Y (0) i+1 = i k=i-2   i l=i-2;l =k s i+1 -s l s k -s l   Y k (40)
Then comes the correction step. Y i+1 needs to be a zero of the function H. The Newton algorithm is used to move

Y (0)
i+1 closer to a zero of H. As Ω is considered as a variable, a path of research needs to be defined. In the pseudo-arclength continuation method, the direction of research

Y i+1 -Y (0) i+1 is orthogonal to the direction of prediction Y (0)
i+1 -Y i , as shown in Figure 10. The (j + 1) th iteration of the correction algorithm is:

J Y H(Y (j) i+1 ) Y (0) i+1 -Y i Y (j+1) i+1 -Y (j) i+1 = -H(Y (j) i+1 ) 0 ( 41
)
The calculation of the jacobian J Y H is made by finite difference where the two values of H are evaluated using the AFT-method.

Once the Fourier coefficients X calculated, the n th harmonic of w is given by the formula

H n = ||A n || 2 + ||B n || 2
, so that the frequency response can be plotted for each harmonic.

Time-integration -Shooting Method

The last simulation method employed in this study is the one presented by Peeters and al. [START_REF] Peeters | Nonlinear normal modes, Part 2: Toward a pratical computation using numerical continuation techniques[END_REF]. An alternative shooting technique coupled with Newmark time marching and an arc length continuation algorithm can also be found in [START_REF] Ibrahim | Modified shooting approach to the non-linear periodic forced response of isotropic/composite curved beams[END_REF]. Problem (32) is written in the following form:

Ż = g(Z) with Z = [W, Ẇ] and g(Z) = [ Ẇ, -D Ẇ -KW + F(W) + F excit ] (42)
The solution of System (42) with the initial conditions Z(t = 0) = z 0 is written as Z(t, z 0 ) in order to exhibit the dependence on the initial conditions. Such a solution is stationary if it is periodic of minimal period T = 1 Ω . A "shooting method" is applied to find out an initial condition z 0 that makes Z(t, z 0 ) to be periodic. Z(t, z 0 ) is calculated by a time-integration method on a discrete time basis of the interval [0, T ]. A fourth order Runge-Kutta integration scheme with an adaptative calculation step control is used. A residual function H is introduced: The frequency response is then the branch of points (z 0 , Ω) satisfying H(z 0 , Ω) = 0. A pseudoarclength continuation scheme is introduced to compute this continuous branch. As presented in Section 5.1, an arclength s is introduced:

H(z 0 , Ω) = Z(T, z 0 ) -z 0 (43) Y i-2 Y i-1 Y i Y (0) i+1 Y i+1 Ω X
∆s i+1 = z 0 (i+1) -z 0 (i) 2 + α Ω (i+1) -Ω (i) 2 (44)
In reference [START_REF] Peeters | Nonlinear normal modes, Part 2: Toward a pratical computation using numerical continuation techniques[END_REF] a tangent predictor is introduced, but a polynomial Lagrangian predictor may be used analogously to Section 5.1. Replacing Y by (z 0 , Ω), the prediction step (40) and the correction steps (41) can be used. The harmonic amplitudes H n can then be computed using a Discrete Fourier Transform algorithm.

Results

In this section, experiments and numerical simulations with ideal or real experimental input are compared.

Simulation with ideal input

First of all, the numerical simulations are compared with experimental data in Figures 11 and12. The results from the two numerical simulation methods are compared for various harmonics. The harmonic curves overlap and differences between each computed harmonic are less than 0.01%. This excellent agreement between time-integration and harmonic balance substantiate the truncature in the Fourier series (p=5) in the harmonic balance method. The computing times of the two methods are also very close and are mainly driven by the continuation algorithm parameters. These simulations are in excellent agreement with experimental data. Figure 11 highlights the differences between the method of multiple scales and numerical methods. The perturbation method is very accurate in the vicinity of F 0 but the accuracy decreases as the frequency distance increases. Focusing on the resonance peak, the analytically calculated peak is higher than the numerically estimated one. The numerical methods give a very accurate solution of the modeled problem thanks to the optimization process. At the resonance peak, the distance between F max and F 0 makes the perturbation method become less accurate than the numerical ones. However, in Figure 11, the analytical results suit experiment better at the resonance peak and could indicate that the perturbation method is the most accurate. This is an artifact due to the updating method. The height of the simulated peak is related to the dissipation parameter µ. This parameter is updated using the perturbation method (see Section 4) and is thus estimated with a bias because of the distance between F max and F 0 . This bias is revealed when the updated model is treated with numerical methods, the numerical results does not suit experiment perfectly. µ is slightly overestimated by the updating method proposed. Figure 12 presents the comparison of measurements with numerical simulation of the harmonic response for various excitation levels. It is clearly seen that the hardening effect is well reproduced via the numerical simulations. These results validate the modeling of the beam with non-ideal boundary conditions and the global strategy proposed in this study, that is based on the coupling between an updating process (using the method of multiple scales) and numerical simulations. Excitation Frequency (Hz) Amplitude of the longitudinal displacement (m) Figure 14: Longitudinal displacement at the clamped ends. u(0) is measured by 1 2 (P2X-P1X) (red crosses). This experimental data is compared with the model simulation (blue plain curve). that a correct estimation of k rot and k bound is necessary to simulate the real structure response. It can be observed that k rot drives the modal participation factor while k bound drives the nonlinear curvature of the frequency response.

The longitudinal displacement at the clamped ends u(0) is also measured during the experiment in order to justify the introduction of the spring k bound . u(0) is computed according to Equation (45) and then compared with experimental data. The experimental signal 1 2 (P2X-P1X) measures u(0). As predicted by Equation (45), this signal's frequency is twice the excitation frequency and its amplitude is compared with simulations in Figure 14. The good agreement illustrated by this figure substantiates the introduction of k bound : the clamped ends shall not be modeled by the ideal condition u(0) = 0.

u(0) = w 2 2 Y 1 (l) 2 1 + lk bound EA -1 l 0 dY 1 dx 2 dx (45)

Simulation with experimental input

In the Harmonic Balance Method presented in Section 5.1, the excitation is assumed to be ideal (F excit = {Γ i a excit } i∈M cos(Ωt)) so that it only contributes to C 1,excit . Anyway, with the same method, the response to a multi-harmonic excitation can be simulated by using coefficients C k,excit = 0. The excitation can thus be the measured input signal, whose harmonic amplitudes are extracted using the method presented in Section 2.3. Then the multi-harmonic response of the model to a multi-harmonic input excitation is computed using Equation (46).

C k,excit = Γ i H k input i∈M and S k,excit = 0 ∀k ∈ [1, p] (46)
Note that in this work only the amplitude of the input signal is used but through S k,excit , the phase information of the input signal can also be taken into account in the simulation. From the authors' experience, changes in the input signal phases do not affect significantly the harmonic amplitudes of the output. The input signal is not independent from the structure response and hysteresis phenomena happen as indicated in Figure 5. The input signal is not the same during the increasing and decreasing swept experiments. For each point of simulation, there should be a choice of experimental input signal (either increasing or decreasing experiment). In this work, simulations reproduce the increasing and decreasing experiments until the jump phenomenon happens. That means that for the first calculation, the continuation method starts from the minimum of the experimental frequency range and increases in frequency until the jump phenomenon. This calculation takes the increasing swept sine experimental data as an input. Analogously, the second calculation starts from the maximum of the frequency range, decreasing in frequency until the second jump phenomenon. This second calculation takes the decreasing swept sine experimental data as an input. The frequency response obtained is not continuous. The continuation method with an ideal input enables to calculate a continuous frequency response beyond singular points but the experimental input cannot be available in the branch beyond singular points.

At the end, the simulation with experimental input can be compared to a simulation with ideal mono-harmonic excitation. These two numerical multi-harmonic responses are then compared with the experimental output, whose harmonic amplitude has also been extracted using the method presented in Section 2.3. These three processes are summed up in Figure 15.

The comparison between experimental data, ideal input simulation and experimental multi-harmonic input simulation is displayed in Figure 16 for the second and third harmonic responses at three levels of excitation. Simulations with the two distinct inputs help understanding the origins of each harmonic response. For the reader's comprehension, it can be noted that the first harmonic experimental response corresponds to the system frequency response. It is predicted by the ideal input simulation, as shown in Figure 12. In this study, the higher order harmonics in the input signal do not affect the first harmonic response. On the contrary, the ideal input simulation does not predict any second harmonic signal whereas it is measured. The experimental second harmonic signal must thus be due to the non-ideal input signal, as previously seen in Section 2.3. As indicated in Figure 5 (a, c and e), the simulations with experimental input confirm this hypothesis: the levels of response predicted are in good agreement with the experimental data for the three levels of excitation. The third harmonic experimental response is explained both by the nonlinearity of the system and by the non-ideal input signal. The simulation with ideal input predicts the third harmonic response in the vicinity of the resonance (see Figure 16, (b, d and f), around 110 Hz). Away from this resonance, the third harmonic experimental response must be due to the non-ideal input. The very good agreement between measurements and multi-harmonic input simulation results substantiates this diagnosis (see Figure 16, (b, d and f), in the ranges [80-100] Hz and [120-160] Hz).

Conclusion

This clamped-clamped beam experiment is made to exhibit, model and simulate nonlinear behaviors of the vibration response of a simple structure with non-ideal boundary conditions. Comparison between experimental and simulation results reveal that the non-ideal boundary conditions are necessary for the modeling. These hypothesis that increase the accuracy of the simulation results are substantiated by complementary measurements that exhibit the predicted longitudinal movement at the boundaries. In this study, three simulation methods adapted to nonlinear vibrations are compared. The method of multiple scales enables to find an analytic equation of the frequency response and thus enables both an easy updating of the model and a good understanding of the influences of each physical parameter on the frequency response. This perturbation method is only accurate in the vicinity of the resonance. Compared to other existing methods, this method could not be easily applied to a more complex structure, such as a finite element model of an industrial assembly. The two numerical methods tested (i.e. Harmonic Balance Method and time-integration method with shooting algorithm) are adapted to any discrete model. The results of the two methods are very similar and in very good agreement with experiments, whatever the frequency and amplitude of excitation. Very similar computational times are found.

Beyond the analysis of the harmonic frequency response, a multi-harmonic comparison of measurements with simulation is performed. A signal processing tool is developed to extract the experimental multi-harmonic frequency response of a structure from a swept-sine experiment. The Harmonic Balance Method enables to compute the multi-harmonic frequency response of a nonlinear model of the structure. The input signal for this simulation can be either an ideal sinusoidal signal or the experimental multi-harmonic input. The comparison between measurements and simulation with these two different inputs enables to distinguish in the frequency response what is due to the system nonlinear response and what is due to the non-ideal input signal.

This study thus outlines the accuracy of the Harmonic Balance Method to predict multi-harmonic responses, and presents a way to interpret the origins of multi-harmonic signals that are often unexpected in the structural frequency responses.

Figure 1 :

 1 Figure 1: A schematic (a) and a picture (b) of the experimental setup which includes 3 accelerometers (P1,P2,A1) and a strain gage (S1).

Figure 3 :

 3 Figure 3: Harmonics' extraction algorithm scheme.

Figure 4 :

 4 Figure 4: Frequency estimation block scheme.

Figure 5 :

 5 Figure 5: Input (a) and output (b) experimental harmonics for an increasing-decreasing swept sine experiment with a constant excitation level of 6 m.s -2 . The arrows indicate the direction of sweep.

Table 1 :Figure 6 :

 16 Figure 6: Schematic of the modeling.

Figure 7 :

 7 Figure 7: First 2 modal shapes.

Figure 8 : 1

 81 Figure 8: Parametric frequency response curve calculated by the method of multiple scales. w low is an arbitrary chosen low amplitude. w low must be low enough to have B w 2 low << 1.

Figure 9 :

 9 Figure9: Comparison of strain gage experimental data (circles) with strain (δε) induced by an updated constant loading (thick straight lines). The comparison is made for the three increasing-decreasing swept sine experiments with excitation levels of 2 m.s -2 (blue), 4 m.s -2 (green) and 6 m.s -2 (red). The offset for the experimental strain values is chosen to set the mean value of the first experiment to zero. Resonance phenomena appear at 100s (increasing sweep) and 480s (decreasing sweep).

Figure 10 :

 10 Figure 10: Schematic of the continuation process. The black curve is the solution. The dashed curve leads to the polynomial prediction. The arrow represents the correction process.

Figure 11 :

 11 Figure11: Comparison of experimental data (red crosses) with simulation by the method of multiple scales (blue dashed curve) and numerical simulations (green curve). The excitation level is 6 m.s -2 . The results of both numerical simulation methods overlap and thus only one curve is plotted.

Figure 12 :

 12 Figure 12: Experimental and simulated harmonic responses at 3 excitation levels. Experimental signal is plotted with crosses. The dashed curves are the simulation results with an ideal mono-harmonic input. The excitation levels are 2 (blue), 4 (green) and 6 (red) m.s -2 .

Figure 13 :

 13 Figure 13: Influence of parameters k bound (a) and k rot (b). Experimental data (red crosses) is compared with numerical simulations. The simulation with updated parameters is plotted on both figures (red thin curve). The blue thick curves represent the asymptotes k bound → ∞ (a) and k rot → ∞ (b). The green dashed curves represent (a) simulation with k bound = 0 and (b) simulation with k rot = 1 2 k rot, updated .

Figure 16 :

 16 Figure 16: Second harmonic responses (a,c,e) and third harmonic responses (b,d,f). The results are obtained with excitation levels of 2 (a,b), 4 (c,d) and 6 m.s -2 (d,e). Experimental signal is plotted with crosses. The dashed curves are the simulation results with an ideal mono-harmonic input. The plain curves are the simulation results with the multi-harmonic experimental input. For the second harmonic, the ideal input simulation results are not plotted because they are < 10 -19 m.

Acknowledgments

The authors gratefully acknowledge the support of the CEA/CESTA teams that helped with this study. They acknowledge in particular the experimental service who provided the measures presented in this paper and also the 2 students Zoé Sigrist and Rémi Andrieux whose works on signal processing enabled to choose and compute the optimization method CMAES.