V Jaumouillé

J-J Sinou

B Petitjean

An Adaptive Harmonic Balance Method for predicting the non linear dynamic responses of mechanical systems -Application to bolted structures

Keywords: Harmonic balance method, bolted joint, nonlinear model, continuation, condensation

Aeronautical structures are commonly assembled with bolted joints in which friction phenomena, in combination with slapping in the joint, provide damping on the dynamic behaviour. Some models, mostly non linear, have consequently been developed and the harmonic balance method (HBM) is adapted to compute non linear response functions in the frequency domain. The basic idea is to develop the response as Fourier series and to solve equations linking Fourier coefficients. One specific HBM feature is that response accuracy improves as the number of harmonics increases, at the expense of larger computational time. Thus this paper presents an original adaptive HBM which adjusts the number of retained harmonics for a given precision and for each frequency value. The new proposed algorithm is based on the observation of the relative variation of an approximate strain energy for two consecutive numbers of harmonics. The developed criterion takes the advantage of being calculated from Fourier coefficients avoiding time integration and is also expressed in a condensation case. However, the convergence of the strain energy has to be smooth on tested harmonics and this constitutes a limitation of the method. Condensation and continuation methods are used to accelerate calculation. An application case is selected to illustrate the efficiency of the method and is composed of an asymmetrical two cantilever beam system linked by a bolted joint represented by a nonlinear LuGre model. The practice of adaptive HBM shows that, for a given value of the criterion, the number of harmonics increases on resonances indicating that non linear effects are predominant. For each frequency value, convergence of approximate strain energy is observed. Emergence of third and fifth harmonics is noticed near resonances both on vibratory responses and on approximate strain energy. Parametric studies are carried out by varying the excitation force amplitude and the threshold value of the adaptive algorithm. Maximal amplitudes of vibration and frequency response functions are plotted for three different points of the structure. Non linear effects become more predominant for higher force amplitudes and consequently the number of retained harmonics is increased. 1

Introduction

The dynamics of mechanical structures is strongly influenced by the presence of riveted or bolted joints in the structure. Indeed structural joints generate energy dissipation through the complex relative motion between two contacting surfaces, commonly referred as frictional slip. Additionally for higher level of excitation, slapping may be encountered. Frictional slip may be analyzed by considering an interface behaviour divided in two cases: micro-slip where part of the interface is slipping; and macro-slip where all the interface slips. Then the frictional energy dissipation observed in the slip zone is responsible for the vibration damping attributed to joints [START_REF] Segalman | Modelling joint friction in structural dynamics[END_REF]. Gaul et al. [START_REF] Gaul | The role of friction in mechanical joints[END_REF] showed that this damping may be larger than material damping and Beards [START_REF] Beards | Damping in structural joints[END_REF] mentioned that up to 90% of the total system damping might be provided by the joints. Thorough reviews about damping in joints may be found in the works of Ungar [START_REF] Ungar | Status of engineering knowledge concerning damping of built-up structures[END_REF], Gaul et al. [START_REF] Gaul | The role of friction in mechanical joints[END_REF] and more recently Ibrahim et al. [START_REF] Ibrahim | Uncertainties and dynamic problems of bolted joints and other fasteners[END_REF].

A better prediction of this damping effect is now an important objective for many aeronautical companies and various complex industrial structures incorporating bolted joints have been investigated [START_REF] Ahmadian | Modelling and updating of large surface-to-surface joints in the awe-mace structure[END_REF][START_REF] Crocombe | Estimating the energy dissipated in a bolted spacecraft at resonance[END_REF][START_REF] Wang | Energy dissipation in spacecraft structures incorporating bolted joints operating in macroslip[END_REF][START_REF] Caignot | Virtual testing for the prediction of damping in joints[END_REF]. Crocombe et al. [START_REF] Crocombe | Estimating the energy dissipated in a bolted spacecraft at resonance[END_REF] established a relationship between energy dissipated in a joint and the transverse excitation force using a 3D FE model of a bolted joint and then used this relationship in conjunction with the simulation of a FE model of a satellite to estimate the energy dissipated in the joints. In the work of Caignot [START_REF] Caignot | Virtual testing for the prediction of damping in joints[END_REF], a micro scale model of bolted joints quantifies in a first step the joint dissipation and an equivalent modal damping is deduced in a second step to perform dynamic analysis of the whole studied structure.

These approaches perform a complex contact analysis giving an insight into the distribution and amount of friction on the interfaces but neglect non linear effects in the global dynamic behaviour of the assembled bolted structure. They need detailed models, often impractical for dynamic analyses of large structures. Hence constitutive models which use a number of degrees of freedom adapted to structural dynamics may be a suitable and computationally efficient alternative. These models can be divided into lumped models and thin layer element theories [START_REF] Ahmadian | Generic element formulation for modelling bolted lap joints[END_REF]. In lumped models, the effect of joint is considered to be concentrated at a single point and the joint model has no dimension. Several models have been proposed: the Valanis model [START_REF] Gaul | Nonlinear dynamics of structures assembled by bolted joints[END_REF], the elasto-slip model [START_REF] Gaul | The role of friction in mechanical joints[END_REF], the LuGre model [START_REF] Gaul | The role of friction in mechanical joints[END_REF], the Iwan model [START_REF] Segalman | Modelling joint friction in structural dynamics[END_REF], the Bouc-Wen model [START_REF] Oldfield | Simplified models of bolted joints under harmonic loading[END_REF], models with Jenkins elements [START_REF] Oldfield | Simplified models of bolted joints under harmonic loading[END_REF] and models integrating a cubic stiffness [START_REF] Ahmadian | Identification of bolted lap joints parameters in assembled structures[END_REF]. The second category, based on thin layer elements, is represented as an element with physical dimensions and specific force-displacement relation. Ahmadian et al. [START_REF] Ahmadian | Generic element formulation for modelling bolted lap joints[END_REF] developed a generic joint element based on a thin layer element approach and Song et al. [START_REF] Song | Simulation of dynamics of beam structures with bolted joints using adjusted iwan beam elements[END_REF] developed an adjusted Iwan beam element incorporating an Iwan model to simulate the dynamics of beam structures.

Most of these models are non linear and require specific methods to compute non linear frequency response functions. In order to compute responses to forced excitation, one of the first methods is time integration. Oldfield et al. [START_REF] Oldfield | Simplified models of bolted joints under harmonic loading[END_REF] applied time integration on a bolted structure to simulate hysteresis loops using a Jenkins element model and a Bouc-Wen model. Other applications on a two beam system were encountered in the works of Gaul et al. [START_REF] Gaul | Nonlinear dynamics of structures assembled by bolted joints[END_REF] and Miller et al. [START_REF] Miller | A two-sided interface model for dissipation in structural systems with frictional joints[END_REF]. Time integration may be inefficient on lightly damped structures because the transient response may take hundreds of forcing periods at the expense of calculation time and disc storage size. Other alternatives like perturbation methods and the Krylov and Bogoliubov method remain limited to a few degrees of freedom.

Heller et al. [START_REF] Heller | Experimental identification of nonlinear dynamic properties of built-up structures[END_REF] applied the Krylov-Bogoliubov method on a non linear system in order to determine equivalent modal parameters and not to compute periodic responses.

In the frequency domain, the harmonic balance method (HBM) is able to compute periodic responses of non linear systems. The basics are to develop the unknown response as a truncated Fourier series and to solve equations linking Fourier coefficients. First mechanical applications can be encountered in the works of Pierre et al. [START_REF] Pierre | Multi-harmonic analysis of dry friction damped systems using an incremental harmonic-balance method[END_REF] on a single degree of freedom dry friction damped system and Ferri et al. [START_REF] Ferri | Frequency-domain solutions to multi-degree-of-freedom, dry friction damped systems[END_REF] on a beam incorporating dry friction. Then a further development of the HBM, named Alternating Frequency Time Domain Method [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF], numerically evaluates the Fourier transform of local nonlinearities of the model and does not require to analytically describe non linear terms. More recently, other approaches have been proposed, notably the Constrained Harmonic Balance Method [START_REF] Coudeyras | A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: The constrained harmonic balance method, with application to disc brake squeal[END_REF] which computes solutions for periodic autonomous systems. For dynamic analyses of bolted joints, Gaul et al. [START_REF] Gaul | The role of friction in mechanical joints[END_REF] used harmonic balance method for the calculation of an equivalent stiffness and viscous damping in an elasto-slip model. Ren et al. [START_REF] Ren | Identification of properties of nonlinear joints using dynamic test data[END_REF] proposed a general technique for identifying the dynamic properties of nonlinear joints using dynamic test data and used multi-harmonic balance method to identify parameters for a friction joint. These two developments compute hysteresis loops and not periodic responses. The work of Ahmadian et al. [START_REF] Ahmadian | Generic element formulation for modelling bolted lap joints[END_REF] developed a non linear generic element formulation for bolted joints and used a cubic non linear stiffness to represent softening non linear effects. Then frequency response curves of the system are calculated with the HBM allowing to include these curves in a minimization procedure in order to identify parameters of the joint. Only prime harmonics were considered due to experimental considerations.

One specific HBM feature is that response accuracy improves as the number of harmonics in the truncated Fourier series increases, at the expense of larger computational time. Therefore only harmonics which lead to a significant contribution on dynamic response must be taken into account for a given precision, and their number can strongly vary on a frequency interval. This key point has been highlighted for bolted joint dynamics by Ouyang et al. [START_REF] Ouyang | Experimental and theoretical studies of a bolted joint excited by a torsional dynamic load[END_REF] who studied an experimental two beam bolted system excited at resonance. By increasing importance of friction in the joints (through an increase of the excitation amplitude), measured hysteresis loops became distorted and superharmonics appeared in the frequency spectra of the responses, showing the importance of considering higher order terms in the Fourier development of the response. Only odd harmonics were present suggesting the possibility to use a cubic stiffness in the bolted joint model [START_REF] Ahmadian | Generic element formulation for modelling bolted lap joints[END_REF][START_REF] Ahmadian | Identification of bolted lap joints parameters in assembled structures[END_REF] and to only consider odd harmonics in the harmonic balance method, usual practice for dry friction system [START_REF] Pierre | Multi-harmonic analysis of dry friction damped systems using an incremental harmonic-balance method[END_REF].

Even though, up to now, no theoretical tool can determine which harmonics are predominant for a non linear system. The present study pursues this investigation by developing a criterion allowing to limit the number of retained harmonics. An approximate strain energy with Fourier coefficients is calculated and its saturation is monitored. This new criterion, based on Fourier coefficients, does not require time integration and may be easily estimated. In order to illustrate the efficiency of the method on a non linear mechanical system, an asymmetrical two cantilever beam system linked by a bolted joint is modelled as application case. The joint model was inspired by the Adjusted Iwan Beam Element (AIBE) developed by Song [START_REF] Song | Simulation of dynamics of beam structures with bolted joints using adjusted iwan beam elements[END_REF]. However, a LuGre model was preferred to an Iwan model present in the work of Song for implementation simplicity. Moreover, formulation of HBM has consequently been adapted to integrate LuGre model internal variables. Analysis of frequency response functions and detailed monitoring of criterion evolution help to assess the validity of this ap-proach. In order to accelerate calculation, a condensation procedure on non linear degrees of freedom is performed by reformulating the HBM equations. Furthermore, criterion has been expressed in this case and compared with case without condensation. This paper is divided into three main sections. The first one deals with the HBM formulation, details condensation and criterion expression. Secondly, the studied system is presented and HBM adaptation to LuGre model is detailed. Finally, result analyses highlight the effect of the harmonic selection process on frequency response functions and on the number of retained harmonics, and a parametric study on the influence of the excitation force is discussed.

HBM Formulation

General Formulation

We consider a discrete mechanical system with nddl degrees of freedoms (dofs) described with its nddl × nddl mass matrix M, stiffness matrix K and damping matrix D. An external periodic force F L (Ω, t) is applied to the system with an angular frequency Ω. System non linearities are considered as a non linear force F N L (X, Ẋ, Ω, t) which depends on degrees of freedom displacements X, velocities Ẋ, angular frequency Ω and time t. The global force F (X, Ẋ, Ω, t) applied on the system may be divided in two parts, the linear external force F L (Ω, t) and the non linear force F N L (X, Ẋ, Ω, t). The governing equation of motion may be written as:

M Ẍ + D Ẋ + KX = F (X, Ẋ, Ω, t) = F L (Ω, t) + F N L (X, Ẋ, Ω, t) (1)
First, we assume a periodic response X(t), which allows to develop the solution as a Fourier series. This development is theoretically infinite so a truncation in the following form is needed:

X(t) = B 0 + m k=1 A k sin( k ν Ωt) + B k cos( k ν Ωt) X(t) = I sin( Ω ν t)I cos( Ω ν t)I . . . sin( k ν Ωt)I cos( k ν Ωt)I . . . [B 0 A 1 B 1 . . . A k B k . . .] T X(t) = T(t)Z (2)
where I is the nddl ×nddl identity matrix,

Z = [B 0 A 1 B 1 . . . A k B k .
. .] T is the (2m+ 1)nddl ×1 vector containing Fourier coefficients, m is the number of harmonics retained for the truncation, ν is an integer used to represent possible subharmonics, and

T(t) = [I sin( Ω ν t)I cos( Ω ν t)I . . . sin( k ν Ωt)I cos( k ν Ωt)I . . .]
is the nddl × (2m + 1)nddl matrix containing trigonometric functions.

The same work is then accomplished for the global force F :

F (X, Ẋ, Ω, t) = C 0 + m k=1 S k sin( k ν Ωt) + C k cos( k ν Ωt) F (X, Ẋ, Ω, t) = T(t) [C 0 S 1 C 1 . . . S k C k . . .] T F (X, Ẋ, Ω, t) = T(t)b (3)
In order to compute velocities and accelerations, we define a frequential derivative operator:

∇ = diag(0 nddl×nddl , ∇ 1 , . . . , ∇ m ) with ∇ k = k ν Ω 0 -I I 0 (4)
Thus we may write:

Ẋ(t) = T(t)∇Z Ẍ(t) = T(t)∇ 2 Z (5)
By replacing Eqn. (2) and Eqn. (5) into Eqn. (1), one obtains:

MT(t)∇ 2 Z + DT(t)∇Z + KT(t)Z = T(t)b (6)
Considering that for a nddl × nddl matrix W and a (2m + 1)nddl × 1 vector Y :

WT(t)Y = T(t)N W Y (7)
with N W = diag(W, W, . . .) (2m + 1)nddl

×(2m + 1)nddl .
Equation ( 6) becomes:

T(t)N M ∇ 2 Z + T(t)N D ∇Z + T(t)N K Z = T(t)b T(t) N M ∇ 2 + N D ∇ + N K Z = T(t)b (8)
Time dependency may be suppressed and a frequency algebraic equation linking Fourier coefficient may be obtained using a Galerkin method which is a projection of the equation on trigonometric functions. Indeed these trigonometric functions define a scalar product:

< f, g >= 1 T T 0 f (t)g(t)dt (9)
Thus we may write:

1 T T 0 T T(t)T(t)dt = 1 2       2I 0 I I 0 . . .       = L (2m + 1)nddl ×(2m + 1)nddl (10)
Applying this scalar product on Eqn. (8) leads to:

1 T T 0 T TT N M ∇ 2 + N D ∇ + N K Z dt = 1 T T 0 T TT b dt L N M ∇ 2 + N D ∇ + N K Z = L b (11)
L is a diagonal matrix so Eqn. [START_REF] Gaul | Nonlinear dynamics of structures assembled by bolted joints[END_REF] may be simplified into a (2m + 1) * nddl equation system:

AZ = b with A = N M ∇ 2 + N D ∇ + N K (12)
A may be expressed in a simpler manner:

A =                   K K -Ω ν 2 M -Ω ν D Ω ν D K -Ω ν 2 M . . . K -k ν Ω 2 M -k ν ΩD k ν ΩD K -k ν Ω 2 M . . .                   ( 13
)
This system is equivalent of finding zeros of a function H : IR (2m+1)×nddl → IR (2 * m+1)×nddl :

H(Z) = A(Ω)Z -b(Z, Ω) (14)
We note that b is dependent on Z and Ω because b corresponds to the Fourier coefficients of F (X, Ẋ, Ω, t).

In the case where no analytical expression may be written between b and Z, an evaluation of the approximate temporal terms X(t) and Ẋ(t) is carried out from an initial value

Z = T [B 0 A 1 B 1 . . . A m B m ]: Z F F T =⇒ X(t) = B 0 + m k=1 (A k sin( k ν Ωt) + B k cos( k ν Ωt)) (15)
It also allows to evaluate temporarily the non linear term F N L (X, Ẋ, Ω, t) and then to deduce Fourier coefficients by a F F T procedure:

F N L (X, Ẋ, Ω, t) F F T =⇒ b N L (Z, Ω) = T [C N L 0 S N L 1 C N L 1 . . . S N L m C N L m ] (16)

Condensation

An additional step can reduce the number of equations to solve. It consists in expressing Fourier coefficients of dofs on which no nonlinearity is applied (called linear dofs) functions of Fourier coefficients of remaining dofs (called non linear dofs) and of Fourier coefficients of linear and non linear forces: First, dofs are reorganized into p linear dofs and q non linear dofs using a boolean transition matrix P:

X = P X p X q = P p P q X p X q (17)
where P p is a nddl×p matrix containing the first p columns of P, P p contains the last q columns of P.

Using the same decomposition for Fourier coefficients,

Z p = [B 0p A 1p B 1p . . . A mp B mp ] T (idem for Z q )
, the following result is obtained:

Z = N Pp N Pq Z p Z q (18)
Note that N Pp and N Pq are respectively (2m + 1)nddl × (2m + 1)p and (2m + 1)nddl × (2m + 1)q matrices. N P and P are both boolean matrices too and T N P N P = I.

Furthermore, the same property is observed for the vector b so that Eqn. ( 12) becomes:

T N Pp AN Pp T N Pp AN Pq T N Pq AN Pp T N Pq AN Pq F pp F pq F qp F qq Z p Z q = b p b q (19)
Z p may be eliminated, and then the system is equivalent of finding zeros of a function

H q : IR (2m+1)×q → IR (2 * m+1)×q : H q (Z q ) = F qq -F qp F -1 pp F pq Z q -b q -F qp F -1 pp b p (20)
If the decomposition is chosen so that no force (linear or non linear) is applied on the p linear dofs, then b p = 0 and:

H q (Z q ) = A q (Ω)Z q -b q (Z q , Ω) (21)
with

A q (Ω) = F qq -F qp F -1 pp F pq .
Finally Fourier coefficients of p linear dofs Z p may be obtained using the relation:

Z p = F -1 pp (b p -F pq Z q ) (22)

Prediction and Correction

For a given frequency Ω, the problem is equivalent to solving a function H(x, Ω) : IR k × IR → IR k with H = H, x = Z and k = nddl or with H = H q , x = Z q and k = q in the condensation case.

When a simulation has to be done on a frequency band [Ω 1 ; Ω 2 ], continuation methods have to be applied to follow the solutions and plot the curve H(x, Ω) = 0. These methods are based on one or more previous points [(x n , Ω n ), (x n-1 , Ω n-1 ), . . .] of the response curve from which a prediction (x

(0) n+1 , Ω (0)
n+1 ) of the next point (x n+1 , Ω n+1 ) is made. Obviously, the closer to the next solution the prediction is, the smaller the number of iterations will be. In order to be able to compute solutions when turning points are present, a curvilinear abscissa s is used. Then, a correction procedure is applied on the prediction, in order to reach after some iterations the next point (x n+1 , Ω n+1 ). In this study no branch points are considered. Then the main prediction and correction methods are presented. The notation y = T [x Ω] ∈ IR k+1 will be used in the following.

Prediction methods -Three prediction methods are presented and illustrated in Fig. 1(a) for y ∈ IR 2 . A given increment ∆s of the curvilinear abscissa is used to calculate the prediction.

Secant Method: the prediction is on the line defined by the two previous points ỹn , ỹn-1 :

y 0 n+1 = y n + ∆s y n -y n-1 y n -y n-1 (23)
Tangent Method: the prediction is on the tangent to the curve at the previous point y n . The direction is given by a unit vector -→ t tangent to the jacobian matrix J y H(y n ) at the point y n so:

y (0) n+1 = y n + ∆s -→ t (24) with J y H(y n ) -→ t = 0 and det J y H(y n ) T -→ t > 0.
Lagrange Polynomial Method: the prediction is on a polynomial P of degree d which reaches the d + 1 previous points [(y n , s n ), . . . , (y n-d , s n-d )]. An analytical description of this polynomial may be easily written by using Lagrange polynomials:

P (s) = n i=n-d          y i n j = n -d j = i s -s j s i -s j          ( 25
)
The prediction is then calculated by evaluating P for the abscissa s n + ∆s:

y (0) n+1 = P (s n + ∆s) (26)
Correction methods -We consider a prediction y

n+1 ∈ IR k+1 of the next solution so the system is not square because the function H offers only k equations. One component of y (often the parameter Ω) has to be fixed or one more equation is added to the system. The correction methods are presented in Fig. 1(b). Correction on x and Ω are noted ∆x and ∆Ω.

Newton Method: one component Ω n+1 is fixed to the prediction value Ω (0) n+1 . Then the system is a square system which can be solved by using a Newton-Raphson procedure. For the i th iteration, the corrected point is:

y (i) n+1 = x(i) n+1 + ∆x Ω n+1 with ∆x = -J x H(x (i) n+1 , Ω n+1 ) H(x (i) n+1 , Ω n+1 ) (27)
Moore-Penrose Method: the vector defined by two consecutive points y

(i+1)
n+1 and y

(i)
n+1 is orthogonal to the kernel of the jacobian matrix J y H(y

(i) n+1 ) of H to the point y (i)
n+1 . The Moore-Penrose matrix inverse defined as W + = T W(W T W) -1 for matrix W is used. The following expression is obtained:

y (i) n+1 = y (0) n+1 + ∆x ∆Ω with ∆x ∆Ω = -J + y H(y (i) n+1 ) H(y (i) n+1 ) (28)
Adaptive step -If a step ∆s used to make a prediction is too large, the number of iterations will be too time-consuming or the solution may not be found. Furthermore, a too large step may lead to difficulties in the vicinity of turning points. So an adaptive step is often appropriate and different methods are available. A widespread tool is to consider the previous number of iterations and to reduce the step when solver takes more iterations than an optimal chosen value. A larger step is chosen when the number of iterations is lower than a minimal value too. Many numerical criteria may be developed to limit the step variation.

(a)

x

y n y (0) n+1 y (0) n+1 Secant Ω y n-1 Tangent x y n Ω y n-2 y (0) n+1 y n-1 Polynomial Prediction (b) x y (0) n+1 y n y n+1 y (i) n+1
Newton Ω

x Ω

y n Moore-Penrose y n+1 y (0) n+1 y (i) n+1

Proposed Criterion

For a given excitation frequency of the system, no theoretical tool exists to determine which harmonics are really predominant. Furthermore the number of necessary harmonics can strongly vary for the studied frequency interval. Some numerical tools have been developed, in particular the work of Laxalde [START_REF] Laxalde | Etude d'amortisseurs non-linéaires appliqués aux roues aubagées et aux systèmes multi-étages[END_REF] who developed a method based on the degree of approximation of the non linearity. This criterion does not take into account the global system behavior and may not be adapted when linear forces are predominant in comparison to non linear forces. Consequently the criterion presented in this section focuses on approximate system strain energy and on its evolution for different numbers of harmonics in the response.

Approximate strain energy -First, for a response developed in a Fourier series, the system strain energy U may be expressed as:

U = 1 2 T X(t)KX(t) = 1 2 T Z T T(t)T(t)N K Z (29)
We may suppress the time dependency by calculating the mean value on one period:

< U >= 1 2 T ZLN K Z (30)
Approximate strain energy for a condensation -When a reduction on q non linear dofs is used, the strain energy expression must be adapted in order to be computed only from Fourier coefficients Z q of the q non linear dofs and to avoid the time consuming step which consists in calculating linear Fourier coefficients Z p .

First some properties must be noted.

For a condensation procedure, the boolean transition matrix P is used. System matrices M (respectively D and K) may be rearranged to correspond to the dofs division [X p X q ] T by using the following

Initialization and Prediction : m = 0 , < Û >= 0, (x conversion M = T PMP (respectively D and K). Similarly to F in Eqn. ( 19), M may be divided in four blocks:

(0) n+1 , Ω (0) n+1
) m = m + 1 HBM Calculation Solution for m harmonics : (x n+1 , Ω n+1 ) m Evaluation of < Ûm >= T xn+1 Kx n+1 Evaluation of criterion ǫ = < Ûm>-< Ûm-1> < Ûm> no Test : ǫ < ǫ threshold yes Solution : (x n+1 , Ω n+1 ) = (x n+1 , Ω n+1 ) m-1 with x =    Z or Z q    and K =    K or Kqq -T Kqp K-1 pp Kpq    .
Mkl = T P k MP l with (k, l) ∈ {p, q} (31)
Furthermore, for a r × s matrix W and a s × t matrix V, N W N V = N WV and T N W = NT W .

Then, if the matrix diag(0 k×k , ∇ 1 , . . . , ∇ m ) with identity matrix as I k×k is named ∇ k , the following formula may be obtained:

∇N P k = N P k ∇ k with k ∈ {p, q} (32)
Finally, the matrix L k is a (2m + 1)k × (2m + 1)k matrix with identity matrix as I k×k and k ∈ {p, q}.

For the computation of criterion in the reduction case, Eqn. (30) becomes:

< U > = 1 2 T ZLN K Z < U > = 1 2 T Z p T Z q T N Pp T N Pq LN K N Pp N Pq Z p Z q < U > = 1 2 T Z p T Z q L p 0 0 L q N PpKPp N PpKPq N PqKPp N PqKPq Z p Z q < U > = 1 2 ( T Z p L p N Kpp Z p + T Z p L p N Kpq Z q + T Z q L q N Kqp Z p + T Z q L q N Kqq Z q ) (33) 10
Then, a more explicit form of matrix F (Eqn. [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF]) must be established:

F kl = T N P k AN P l (34)
Introducing the expression of matrix A of Eqn. ( 12) leads to:

F kl = T N P k N M ∇ 2 + N D ∇ + N K N P l F kl = N Mkl ∇ 2 l + N Dkl ∇ l + N Kkl (35)
with (k, l) ∈ {p, q}. For simplicity, linear and non linear Fourier coefficients are supposed to be linked in a static case (Ω = 0) so that ∇ l = 0 l×l . By using this assumption in Eqn. [START_REF] Ouyang | Experimental and theoretical studies of a bolted joint excited by a torsional dynamic load[END_REF], Z p may be expressed as:

Z p = N -1 Kpp b p -N Kpq Z q (36)
If the decomposition is chosen so that no force (linear or non linear) is applied on the p linear dofs, then b p = 0 and:

Z p = -N -1 Kpp N Kpq Z q (37)
By replacing Eqn. (37) into Eqn. (33):

< U > = 1 2 ( T Z q T N Kpq T N K-1 pp L p N Kpp N K-1 pp N Kpq Z q + T Z q T N Kpq T N K-1 pp L p N Kpq Z q + T Z q L q N Kqp N K-1 pp N Kpq Z q + T Z q L q N Kqq Z q ) < U > = 1 2 T Z q L q N Kqq-T Kqp K-1 pp Kpq Z q (38)
This equation is very similar to Eqn. (30) obtained for the non reduced case. Matrix Kqq -T Kqp K-1 pp Kpq acts as reduced stiffness matrix on the q non linear dofs.

Criterion ǫ -The criterion ǫ developed in this section is computed for a given frequency and the first resolution is performed for one harmonic. Then the relative difference between two consecutive values of strain energy is evaluated. The first value is obtained for m harmonics and the second for m + 1 harmonics. The increase is stopped when ǫ becomes less than a threshold chosen by user. Algorithm is detailed in Fig. 2. As matrices L and L q are diagonal constant block matrices and as algorithm starts for one harmonic, studying strain energy saturation is equivalent to studying saturation of an approximate quantity:

< Û > = T ZN K Z without condensation < Û > = T Z q N Kqq-T Kqp K-1
pp Kpq Z q with condensation (39)

Finally, when the convergence rate of the strain energy turns out to be non smooth, the method may stop before saturation. For example, for dry friction systems, only odd harmonics appear. However, as shown later, this drawback may be avoided removing all even harmonics in the calculation. For other cases, for example when the 1 st and 5 th harmonics responds, this drawback constitutes a limitation of the method.

3 Application Case

Two Beam System and Joint Model

A two cantilever beam system linked by a bolted joint is considered for simulations in order to illustrate the selection process of harmonics and is shown in Fig. 3 For further details on geometric non linearities, the reader can refer to Sze et al. [START_REF] Sze | The incremental harmonic balance method for non linear vibration of axially moving beams[END_REF] who applied HBM on a non linear beam. Bolted joint is represented with a 3.5cm long element described by a mass elementary matrix of a B21 beam element and a non linear stiffness matrix considered as an external force. These external forces are the two moments M 1 and M 2 and the two forces T 1 and T 2 . A Rayleigh damping is calculated by using mass and stiffness matrices of a monolithic beam so that damping has a value of 0.1% for frequencies of 0.24kHz and 1.14kHz. These two frequencies correspond to the second and fifth modes. The value of 0.1% is representative of real structures and remains sufficiently low to keep a significative non linearity impact. Indeed, a material damping increase leads to reduce vibration amplitudes and consequently higher order harmonics amplitudes so adaptive harmonic balance method will retain a smaller number of harmonics. Non linear stiffness model of bolted joint presented in Fig. 3(b) is largely inspired by a previous work of Song [START_REF] Song | Simulation of dynamics of beam structures with bolted joints using adjusted iwan beam elements[END_REF] and uses a model called Adjusted LuGre Beam Element (ALBE) by analogy with Song's work. Frictional slip and slapping constitutes the two main non linear phenomena involved in the joint interface [START_REF] Ahmadian | Generic element formulation for modelling bolted lap joints[END_REF]. However, for the sake of simplicity, slapping which can lead to higher order harmonics has not been considered in this study. Frictional slip is here considered with a non linear model integrating a LuGre model leading to odd harmonics. As discussed in section 2.4, even harmonics have to be removed in the HBM calculation so that smooth convergence rate of strain energy is observed on only odd harmonics. The basic idea is to replace stiffnesses of a linear beam element by a parallel combination of a LuGre model and a residual stiffness k a,i , i ∈ 1, 2 characteristic of a bolted joint [START_REF] Gaul | Nonlinear dynamics of structures assembled by bolted joints[END_REF]. The element has two rotational dofs R1 and R2 and two translational dofs T 1 and T 2. h and L are respectively section height and element length. Spring elongations ∆ 1 and ∆ 2 have to be considered to express relation between non linear force F N L,ALBE and element dofs:

∆ 1 = L 2 (R 1 + R 2 ) + (T 1 + T 2) and ∆ 2 = h 2 (R 1 -R 2 ) (40)
Consequently each LuGre force f LuGre,i , i ∈ 1, 2 depends on ∆ i elongation but also on internal variable value ζ i and its derivative ζi . It may be written as:

f LuGre,i (∆ i , ∆i , ζ i , ζi ) = σ 0i ∆ i + σ 1i ζi + α 2i ∆i (41) ζi = ∆i - σ 0i α 0i + α 1i e -∆i v 0i 2 ∆i ζ i (42)
The combination of one LuGre model and one spring provides a force which takes the following form:

f i (∆ i , ∆i , ζ i , ζi ) = f LuGre,i (∆ i , ∆i , ζ i , ζi ) + k a,i ∆ i (43)
Stiffness decrease during microslip regimes may be represented by using a coefficient γ i ∈ [0; 1] which links the residual stiffness k a,i , the LuGre model stiffness parameter σ 0i and the equivalent linear element stiffness k i . The relative equations are

σ 0i = (1 -γ i )k i et k a,i = γ i k i .
Forces F 1 , F 2 and resulting moments M 1 , M 2 may be expressed as:

     F 1 M1 F 2 M2      =       f 1 (∆ 1 , ∆1 , ζ 1 , ζ1 ) L 2 f 1 (∆ 1 , ∆1 , ζ 1 , ζ1 ) + h 2 f 2 (∆ 2 , ∆2 , ζ 2 , ζ2 ) -f 1 (∆ 1 , ∆1 , ζ 1 , ζ1 ) L 2 f 1 (∆ 1 , ∆1 , ζ 1 , ζ1 ) -h 2 f 2 (∆ 2 , ∆2 , ζ 2 , ζ2 )       (44)
Thus from Eqn. (40) forces F N L,ALBE may be written as a function of ALBE dofs:

F N L,ALBE (      T 1 R 1 T 2 R 2      ,      T 1 R 1 T 2 R 2      , ζ 1 ζ 2 , ˙ ζ 1 ζ 2 ) =      F 1 M1 F 2 M2      (45)
The two equivalent linear element stiffnesses k i are obtained by the following relations: Then the system is excited with a harmonic excitation F L of 42N with a pulsation Ω. Load is applied on the longest beam on translational dof 31. By using formulation introduced in Eqn.

k 1 = 12 EI L 3 = 1, 43.
(1), one may write the governing equation of motion to which two equations have to be added. These additional equations describe LuGre model internal variables ζ 1 , ζ 2 evolution:

M Ẍ + D Ẋ + KX = F L (Ω, t) -F N L,ALBE ( X Ẋ , ζ 1 ζ 2 , ζ1 ζ2 , Ω, t) (47) ζi = ∆i - σ 0i α 0i + α 1i e -∆i v 0i 2 ∆i ζ i f or i ∈ {1, 2} (48)

Adaptation of HBM Formulation to LuGre Model

In the studied case, two equations are added to the equation of motion and integrate non linear terms. Thus an adaptation of HBM formulation becomes necessary. To do so, the two internal variables ζ 1 and ζ 2 are developed as a Fourier series in the same way as X. By inserting the truncature 48), the same Galerkin method is applied on the resulting equation and one obtains:

ζ i (t) = T(t)Z ζ i into Eqn. (
0 = ζi -     ∆i - σ 0i α 0i + α 1i e -∆i v 0i 2 ∆i ζ i     0 = T(t)∇Z ζ i -T(t)b ζ i (Z, Z ζ i , Ω) 0 = N I 1×1 ∇Z ζ i -b ζ i (Z, Z ζ i , Ω) with i ∈ {1, 2} (49)
ζ = [Z ζ 1 Z ζ 2 ].
HH is a function from IR (2m+1)×(nddl+2) to IR (2 * m+1)×(nddl+2) .

HH(Z, Z

ζ ) = H(Z, Z ζ ) = A(Ω)Z -b(Z, Z ζ , Ω) C(Z, Z ζ ) = N I 2×2 ∇Z ζ -b ζ (Z, Z ζ , Ω) (50)
In case of condensation, we may define similarly a function HH q (Z q , Z ζ ) from IR (2m+1)×(q+2) to IR (2 * m+1)×(q+2) :

HH q (Z q , Z ζ ) = H q (Z q , Z ζ ) = A q (Ω)Z q -b q (Z q , Z ζ , Ω) C(Z q , Z ζ ) = N I 2×2 ∇Z ζ -b q,ζ (Z q , Z ζ , Ω) (51)
Finally, calculation of approximate strain energy < Û > stay the same as Eqn.(33) and Eqn. (38) because only Fourier coefficients of physical dofs are used to quantify strain energy.

Results

In the following, shown results have been calculated on the frequency band [0 -2.3] kHz with a curvilinear abscissa and an adaptive step to better describe resonance peaks. The Moore-Penrose method is used for correction at each iteration and prediction is made with Lagrange polynomials of degree 2.

Non linear effects on dynamic responses

The maximal amplitude of vibration obtained for the linear case and for the two different non linear cases (1 harmonic case and adaptive case) are presented in Fig. 4(a,b,c) for the three dofs of the system. The 1 harmonic case refers here to the classical HBM with one harmonic and does not refer to an adaptive algorithm. The first dof 7 is located on the left beam, the second is a translational dof 19 corresponding to the R1 dof of the ALBE model, and the third one is the dof 59 on the right beam.

Figure . 3(a) details the dof position on the system. The threshold value for the relative variation of the approximate strain energy has been fixed to 3% and excitation force has an amplitude of 42N.

The linear case shows seven modes on the studied frequency band. First non linearity effects are significant on the one harmonic response and result in two phenomena: a reduction of resonance peak amplitudes which reflects the damping from the joint and a modal softening which reflects the joint stiffness decrease. Thus some modes have important frequency shifts and even significant distortions, notably the second, fifth and sixth modes. Frequency shifts and vibration amplitude reductions have the same order of magnitude for all the considered dofs. Differences are observed between the 1 harmonic curve and the adaptive algorithm near resonances especially near the fifth and sixth modes where the shape of peaks differs. Far from resonance peaks, adaptive algorithm give the same results as the 1 harmonic calculation.

Figures 4(d,e,f) present a zoom on the 1.1kHz resonance peak. Three non linear cases corresponding to a calculation with 1, 3 and 11 harmonics are compared with the adaptive HBM curve. Near the resonance, the adaptive HBM remains close to the response with 11 harmonics. However, Figure 5, which plots the number of harmonics over the whole frequency band, shows that the number of harmonics reaches only 7 harmonics at most, revealing satisfying convergence of the method. The number of used harmonics may vary from 1 to 11 for resonance peaks but only 1 harmonic is necessary elsewhere. It shows that non linear phenomena are more pronounced on resonances. The maximum number of harmonics has been fixed to 15 in this case.

Contribution of each harmonic (1 st ,3 rd and 5 th harmonic) is shown in Fig. 6(a,b,c) by plotting the Fourier coefficient modulus of the 1 st ,3 rd and 5 th harmonics of the vibration response (||A k B k || for harmonic k and for one dof into Eqn. ( 2)) near the 1.1kHz resonance peak. Analyses are still performed on the three dofs 7, 19, and 59. For the all three dofs, the 1 st harmonic contribution remains larger than the 3 rd and 5 th harmonic contributions. However, analysis of the ratios 3 rd /1 st (Fig 6 (d,e,f)) and 5 th /1 st (Fig 6 (g,h,i)) shows that importance of third anf fifth harmonics increases near resonances and may represent up to 10% for the third harmonic and up to 2% for the fifth harmonic. This observation may be linked with the work of Ouyang et al. [START_REF] Ouyang | Experimental and theoretical studies of a bolted joint excited by a torsional dynamic load[END_REF] who found emergence of third and fifth superharmonics which represented about respectively 2% and 0.7% of the first harmonic. It may also be noted that third and fifth harmonics are less predominant for dof 7 than for dofs 19 and 59.

A similar analysis on the approximate strain energy may be carried out by considering the contribution of the order k as being the term 1 2 T Z k KZ k . Z k refers here to the contribution of the order k to the Fourier coefficient vector Z of the vibration response. Results are plotted in Fig 7(a). Ratios 3 rd /1 st and 5 th /1 st are shown in Fig 7(b,c). The same tendency that for maximal amplitude of vibration with a peak near resonance frequencies is observed. Moreover, ratios have the same order of magnitude with a maximum value of 12.5% for the third harmonic and of 2% for the fifth harmonic, revealing that approximate strain energy behaves like a global indicator of each dof behaviour.

Approximate Strain Energy Saturation

Approximate strain energy < Û > is presented over all the frequency band for three cases on Fig. 8(a). The first case shows the shape of < Û > for a linear case and a peak is observed for each resonance frequency. The other two cases correspond to non linear calculations with one harmonic and with an adaptive number of harmonics. Non linear effects decrease vibration amplitude of the system so that peaks on approximate strain energy are attenuated. As observed on the maximal amplitude of vibration curves, peaks are shifted to the left. Moreover, differences between linear and non linear cases are predominant near resonance frequencies and adaptive algorithm curve differs from one harmonic curve showing an increase in the required number of harmonics. A zoom near the 1.1kHz peak is made on Fig. 8(b) in order to show convergence of the approximate strain energy quantity. Results are presented for 1,3 and 11 harmonics and for the adaptive case. Saturation is observed and adaptive case stays close to 11 harmonic curve even if no more than 9 harmonics are used. Finally, it has to be noted that strain energy saturation is directly tested on odd harmonics due to the presence of dry friction in the model and considered harmonics have a monotonous decrease of their amplitude. It constitutes a limitation of the method which cannot deal with non consecutive predominant harmonics (for example system with 1, 3 and 11 predominant harmonics).

Influence of force excitation

First simulations were carried out with an excitation amplitude of 42N. Influence of excitation force amplitude on non linear effects is now investigated by varying amplitude with the values 6N, 12N, 24N, 42N, 66N. The maximal amplitude of vibration is presented for the three considered dofs in Responses are computed by using the adaptive algorithm. First, we note that an increase in excitation amplitude results in larger vibration amplitudes, and this for all the three dofs. Then, the main notable non linear effect is an increase in modal softening for larger excitation amplitude, as shown by left shifts of resonance peaks. It clearly shows a relationship between modal softening and vibration amplitudes and this dependence is non linear as previously noticed by Ungar [START_REF] Ungar | Status of engineering knowledge concerning damping of built-up structures[END_REF]. For larger amplitudes, this modal softening becomes less remarkable revealing the beginning of macro-slip and so the stabilization of the contact stiffness. In order to compare the five non linear cases, frequency response functions (FRFs) are computed by dividing the maximal amplitude of vibration by the excitation force amplitude for each frequency value. Results, which are very similar for the three considered dofs, are presented in the particular case of the dof 19 in (c) and for the seventh peak (d). First we note that modal softening is still observed on resonance peaks and FRF maximal value decreases when excitation amplitude increases, reaching a minimal value between the 24N and 42N cases before a new increase for larger excitation amplitudes, except for the 2.2kHz peak where only a decrease is observed. Damping provided by joints is non linear and amplitude dependent and this may explain this behaviour. Ouyang et al. [START_REF] Ouyang | Experimental and theoretical studies of a bolted joint excited by a torsional dynamic load[END_REF] observed the same phenomenon on the relation between energy dissipation and excitation amplitude for a torsional joint and notices that, as excitation amplitude becomes larger, micro-slip increases.

Excitation amplitude has an influence on non linear effects and consequently on the number of harmonics used by the adaptive algorithm. Figure 12(a,b,c,d,e) presents the number of obtained harmonics for the five different excitation values. A zoom on the 1.1kHz peak is done. For the lowest amplitude value, the adaptive HBM only requires one harmonic whereas for the highest amplitude it uses nine harmonics. This increase is progressive and the frequency interval for which more than one harmonic are needed widens. Left shifts of resonances are also noticeable because the frequency for which the number of harmonics is maximum decreases as excitation amplitude becomes higher.

Influence of threshold value

Three threshold values of 1%, 3%, 5% have been tested for the adaptive algorithm. Figure 13(a,b,c) shows the evolution of the number of harmonics for the two peaks near 1.1kHz and 1.5kHz and for these three threshold values. Unsurprisingly, an increase of the selected number of harmonics is observed when the threshold value decreases because more harmonics are needed to reach this precision. Threshold value still remains a value to be determined by user.

(a)

Analysis of general and reduced criteria

An expression of approximate strain energy has been established in Eqn. (38) when a condensation procedure is used in order to avoid the calculation of linear Fourier coefficients for each step of the adaptive algorithm. This expression is based on an assumption that linear and non linear Fourier coefficients are linked in a static case allowing to obtain an approximate strain energy described only with the reduced stiffness matrix and not with mass and damping matrices. A calculation with this reduced criteria has been carried out and the evolution of the approximate strain energy for the two condensation and general cases are overlaid in Fig 14(a). The two quantities clearly differ for frequencies higher than that of the first resonance, invalidating the assumption used for the calculation of the strain energy. However Fig. 14(b), which plots the evolution of the number of obtained harmonics when the reduced criterion is used, shows that the result of the adaptive algorithm is very close to the result for the general case presented in Fig. 5. This may be explained by the fact that the used criterion computes a relative quantity and observes strain energy saturation by evaluating a relative difference between two consecutive values. The use of this criterion may represent an alternative for an adaptive HBM calculation.

(a)

Conclusion

This article develops a new adaptive harmonic balance method which selects the number of harmonics at each frequency value for a mechanical system integrating localized non linearities. A two beam system linked by a bolted joint is chosen for application case. HBM formulation is combined with a reduction on non linear dofs of the system and a simulation on a frequency band is carried out by using a prediction method based on Lagrange polynomials and a correction method based on Newton and MoorePenrose methods. In order to adapt the number of harmonics taken into account at each frequency, an adaptive algorithm has been developed. The method computes an approximate strain energy from Fourier coefficients of the response and observes its saturation by evaluating the relative difference between two consecutive cases corresponding to two different number of harmonics. This new criterion, based on Fourier coefficients, does not require time integration and may be easily estimated. In a condensation case, criterion formulation is also expressed. Furthermore geometric non linearity due to axial forces in the clamped-clamped beam are not taken into account and non linear effects in the joint consider only frictional slip. Slapping is not modelled in this study. Slip in the bolted joint element is represented by a LuGre model which leads to adapt the HBM formulation in order to develop internal variables as Fourier series.

Results show that one harmonic is sufficient to give a satisfactory approximation of the response away from resonances and is necessary to highlight non linear effects such as damping of resonance peaks and modal softening. Indeed the dynamic behaviour is strongly modified compared with the linear case. Moreover, adaptive HBM shows that, for a given threshold value of the criterion, the number of harmonics may increase on resonances indicating that non linear effects are predominant. The evolution of the approximate strain energy shows that a peak is observed near each resonance and saturation of this quantity is noted when the number of harmonics increases. However, calculation is performed only on odd harmonics due to dry friction leading to a smooth convergence rate of the strain energy on tested harmonics. This condition constitutes a limitation of the method which cannot deal with non consecutive predominant harmonics (for example system with predominant harmonics 1, 3 and 11). Analysis of each harmonic contribution notices the emergence of third and fifth harmonics both on the response and on approximate strain energy near resonances, showing the global characteristic of the criterion based on approximate strain energy. In order to obtain a wider range of harmonics and to model a more physical bolted system, slapping and geometric non linearities could be considered for further work. A coherent behaviour is noticed when threshold value varies because more harmonics are needed to reach the given precision when threshold value of the adaptive algorithm is decreased. A parametric study is carried out by varying the excitation force amplitude. Vibration amplitude increases with higher force amplitude because non linear effects, notably micro slip in the joint, become more pronounced. Modal softening and damping depends on vibration amplitude and this dependency is non linear. Maximum of frequency response functions for each resonance depends non linearly on excitation amplitude and may reach a minimum value for an intermediate excitation amplitude. The number of needed harmonics becomes larger for increasing amplitudes underlining the predominance of non linear effects.

Figure 1 :

 1 Figure 1: Prediction and Correction Methods.

Figure 2 :

 2 Figure 2: Algorithm for criterion ǫ.

Figure 3 :

 3 Figure 3: Application Case (a) Two beam system and (b) bolted joint model.

  10 6 N/mm et k 2 = 4 EI Lh 2 = 8, 92.10 5 N/mm (46) Other parameters are deduced by analogy with Shiryayev's work [25], namely γ 1 = γ 2 = 0.1078,α 01 = α 02 = 81.9N, σ 11 = σ 12 = α 11 = α 12 = α 21 = α 22 = 0. Finally σ 01 = 1, 27.10 6 N/mm, σ 02 = 7, 96.10 5 N/mm, k a,1 = 1, 55.10 5 N/mm, k a,2 = 9.62.10 4 N/mm.

  where b ζ i (Z, Z ζ i , Ω) are Fourier coefficients of the non linear term. Consequently the problem is equivalent to finding zeros of a function HH(Z, Z ζ ) depending on Fourier coefficients of X and ζ = [ζ 1 ζ 2 ]. These coefficients are named Z and Z

Fig 9 .

 9 Zooms for the dof 19 are done for the first three peaks in Fig 9(a), for the fifth peak in Fig 9(b) and for the seventh peak in Fig 9(c).

  Fig 11 for all the frequency band (a), for the first three peaks (b), for the fifth peak

Figure 4 :

 4 Figure 4: Maximal amplitude of vibration for dofs 7(a,d), 19(d,e), 59(c,f) and for different linear and non linear cases: (a-c) over all the frequency band; --(linear, HBM 1 harm.); -(non linear, HBM 1 harm.); ... (non linear, adaptive HBM) (d-f) zoom on 1.1kHz resonance; -(non linear, HBM 1 harm.); -.-(non linear, HBM 3 harm.); ... (non linear, adaptive HBM); -.-(non linear, HBM 11harm.)

Figure 5 :

 5 Figure 5: Evolution of the number of harmonics: -number of harmonics; ... maximal amplitude of vibration for dof 19.

Figure 6 :

 6 Figure 6: Contribution of 1 st ,3 rd and 5 th harmonics on maximal amplitude of vibration and associated ratios for dofs 7(a,d,g), 19(b,e,h), 59(c,f,i): (a-c) Fourier coefficient modulus of the: -1 st harmonic; --3 rd harmonic; ... 5 th harmonic; (d-f) ratio 3 rd /1 st ; (g-i) ratio 5 th /1 st .

Figure 7 :Figure 8 :

 78 Figure 7: Contribution of 1 st ,3 rd and 5 th harmonics on approximate strain energy and associated ratios: (a) contribution of the: -1 st harmonic; --3 rd harmonic; ... 5 th harmonic; (b) ratio 3 rd /1 st ; (c) ratio 5 th /1 st .

Figure 9 :

 9 Figure 9: Influence of excitation force amplitude on maximal amplitude of vibration for dofs 7(a), 19(b), 59(c): ... 6N; ---12N; -.-24N; --42N; -66N.

Figure 10 :

 10 Figure 10: Influence of excitation force amplitude on maximal amplitude of vibration for dof 19: (a) Zoom on 3 first resonances, (b) Zoom on 1.1kHz resonance, (c) Zoom on 2kHz resonance: ... (replaced by +++ on (c)) 6N; ---12N; -.-24N; --42N; -66N.

Figure 11 :

 11 Figure 11: Influence of excitation force amplitude on frequency response functions for dof 19: (a) over all the frequency band, (b) Zoom on 3 first resonances, (c) Zoom on 1.1kHz resonance, (d) Zoom on 2kHz resonance : ...(replaced by +++ on (d)) 6N; ---12N; -.-24N; --42N; -66N.

Figure 12 :

 12 Figure 12: Influence of excitation force amplitude on the number of harmonics, zoom on the 1.1kHz resonance: (a) 6N, (b) 12N, (c) 24N, (d) 42N, (e) 66N.

Figure 13 :

 13 Figure 13: Influence of threshold value on the number of harmonics, zoom on the 1.1kHz and the 1.5kHz resonances: (a) 5%, (b) 3%, (c) 1%.

Figure 14 :

 14 Figure 14: (a) Approximate strain energy: -general expression; --condensation case expression (b)Number of harmonics obtained with the condensation case expression.

  (a). Beams are made of aluminium 7075 Al and the section is rectangular (5.1cm × 2.5cm). The two beams have different lengths (34.7cm and 84.7cm) in order to avoid a symmetric behaviour. Beams are modelled with Abaqus software with two dimensional B21 beam elements which use a Timoshenko formulation. 10 (respectively 25) elements are used for the 34.7cm (respectively 84.7cm) beam. Due to axial forces, geometric non linearities may appear in this clamped-clamped beam system for high level of excitation and may be simultaneously present with joint non linearity. However this model is focused on localized non linearities and study has been limited to joint non linearities.


F Chevillot

J.-J Sinou

N Hardouin

L Jézéquel

Effects of damping on the speed of increase and amplitude of limit cycle for an aircraft braking system subjected to mode-coupling instability

Keywords: friction-induced vibration, damping, limit cycles, nonlinear dynamics

A nonlinear model of an aircraft braking system is presented and used to investigate the effects of damping on the stability in Chevillot et al. (Arch Appl Mech 78(12):949-963, 2008). It has been shown that the addition of damping into the equations of motion does not lead systematically to the stabilization of the system. In the case of a mode-coupling instability, there is indeed an optimal ratio between the modal damping coefficients of the two modes in coalescence, that maximize the stable area. But the stable area is not a sufficient criterion. In dynamics, the amplitude of the vibrations and the transient behavior characterized by the speed of increase of the oscillations are best indicators. In this paper, the same nonlinear model of the aircraft braking system is used to compute time-history responses by integration of the full set of the nonlinear dynamic equations. The aim of the study is to evaluate the effects of damping on the nonlinear dynamics of the brake. It is shown that damping may be very efficient to significantly reduce and slow down the increase of the friction-induced vibrations. But, in the same way as for the stability area, there exists a value of the damping ratio that optimizes the effects of damping.

Introduction

The problem of unstable vibrations in mechanical systems with rotating and sliding parts has received the attention of many investigations, and is still and active field of research in dynamics. Solving potential vibration problems requires the consideration of the stability analysis and the determination of the time-history response. A stability analysis as presented in The destabilization paradox applied to friction-induced vibrations in an aircraft braking system [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] allows determining the characteristics (frequency, real part, mode shape and critical friction coefficient) of the instabilities that may appear during the braking by extraction of the eigenvalues of the linearized system at the equilibrium point. In particular it allows knowing if for a given set of parameters the system is stable or not. But if not, it gives no information on the amplitude of the oscillations or on the nonlinear transient behavior. To achieve this, the complete nonlinear response of the system has to be computed.

In this paper, the model developed in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] is used to investigate the effects on damping on the nonlinear dynamical behavior of an aircraft braking system. The results of [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] indicate that the addition of damping in the equations of motion alter the stable-unstable boundary. In some cases, a destabilizing effect was shown: damping is not always beneficial. The paper also shows that for this particular system there exists an optimal structural damping ratio between the stable and unstable modes for which the stable region is the largest.

The effects of damping on non-linear non-conservative systems and the highlight of a destabilizing paradox was first studied by Ziegler [START_REF] Ziegler | Die stabilitätskrieterien der elastomechanik[END_REF]. The study was focused on a double pendulum under a follower force. He showed that the critical load of the damped system was significantly smaller than that of the undamped system. The fact that damping can be harmful to a system is not intuitive and runs counter to the generally accepted idea of the effects of damping. Many investigators have subsequently looked into this phenomenon. In particular, Bolotin [START_REF] Bolotin | Non-conservative problems of the theory of elastic stability[END_REF] showed that the ratio between the damping coefficients plays a key role in the destabilizing effect. Later on, Herrmann [START_REF] Herrmann | Stability of equilibrium of elastic systems subjected to nonconservative forces[END_REF], Walker [START_REF] Walker | A note on stabilizing damping configurations for linear non-conservative systems[END_REF], O'Reilly et al. [START_REF] O'reilly | Some aspects of destabilization in reversible dynamical systems with application to follower forces[END_REF] and recently Kirillov and Seyranian [START_REF] Kirillov | Stabilization and destabilization of a circulatory system by small velocity-dependent forces[END_REF][START_REF] Kirillov | The effect of small internal and external damping on the stability of distributed non-conservative systems[END_REF] worked on the determination in explicit form of the stabilizing damping configurations for a large class of non-conservative systems. In particular in [START_REF] Kirillov | Stabilization and destabilization of a circulatory system by small velocity-dependent forces[END_REF] necessary and sufficient conditions for the matrices of velocity-dependent forces to stabilize an unperturbed circulatory system have been expressed. Applied to the Ziegler pendulum, it was found that the optimum damping ratio that gives the larger critical load is far from unity. In wood cutting machines, Gallina and Trevisani [START_REF] Gallina | On the stabilizing and destabilizing effects of damping in wood cutting machines[END_REF] demonstrated that similar damping coefficients cannot destabilize the initially stable system. Using a finite element model Fritz [START_REF] Fritz | Effects of damping on brake squeal coalescence patterns -application on a finite element model[END_REF] came to the same conclusion for an automotive brake squeal model. An asymmetry in the damping distribution tends indeed to make the system unstable for a lower value of a the friction coefficient. However, Sinou and Jézéquel [START_REF] Sinou | Mode coupling instability in friction induced vibrations and its dependency on system parameters including damping[END_REF] for example studied the effects of damping with a nonlinear two-degree-of-freedom model with friction and found that an optimal damping distribution not always exists. The damping coefficients, their ratio, and the pulsation ratio have to be considered. Hervé et al. [START_REF] Hervé | Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects[END_REF] used a two-degree-of-freedom nonlinear model of clutch squeal to investigate the influence on circulatory and gyroscopic actions. While for the purely circulatory system proportional damping is found to be the optimum, taking into account the gyroscopic effects makes a lot more complex the effects of damping. Moreover, most of these studies focused on the effects on the stability region. However, because a stability analysis is based on the nature of the eigenvalues of the linearized equations of motion, it is not obvious that the effects on the nonlinear system are the same. Aware of that, recent papers (Kounadis [START_REF] Kounadis | On the paradox of the destabilizing effect of damping in nonconservative systems[END_REF], Bolotin et al. [START_REF] Bolotin | Effect of damping on the postcritical behaviour of autonomous non-conservative systems[END_REF], Sinou and Jézéquel [START_REF] Sinou | Mode coupling instability in friction induced vibrations and its dependency on system parameters including damping[END_REF] for instance) included in their studies the analysis of the steady-state amplitudes. It is motivated by finding that regions considered unstable by a stability analysis turn out to be stable when a non-linear dynamic analysis is performed (Kounadis [5]).

In this paper, the model developed in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] is briefly exposed. Then, the analysis of the stability of the system is presented. While Chevillot et al. [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] deal with three low-frequency mode-coupling instabilities calculated, this study focuses on the second instability, for which the evolution of the stability area versus the friction coefficient and the damping ratio is reminded. Time-history simulations with variations in the friction coefficient and damping configuration are performed to analyze their effects on the transient behavior and the limit cycle amplitude of the instability. The main purpose of the present paper is to find if the conclusions expressed in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] about the stability regions are similar when a nonlinear dynamic analysis is performed. The results show that damping is a key factor very efficient to reduce frictioninduced vibrations, but has to be used cautiously to avoid misuse. In particular, the destabilizing effect and the optimum damping distribution highlighted in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] for the stability analysis are extended to the nonlinear dynamics of the system. Based on experimental observations, an analytical model of an aircraft braking system given in Fig. 1a has been developed and presented in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] in order to reproduce whirl and squeal instabilities. The view of the parametric model using 70 degrees of freedom in Fig. 1b allows seeing the correspondence between the model and the system. The multi-stage brake is represented by a single rotor and stator with an effective brake friction coefficient of µ eq = 2nµ where n is the number of rotors, producing 2n interfaces in contact between the stators and rotors. It is assumed that the friction coefficient is uniform along the rubbing surface and that the rotor and stator friction surfaces are always in contact. Friction is modeled by a Coulomb approach: the tangential force is related to the normal force by the friction coefficient (for more details see [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF]). Experimental tests show that the behavior of the rotor/stator stack in compression is nonlinear. The analytical model uses a third order polynomial to express the load-deflection relationship [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF].

The global nonlinear forces and moments due to friction are then expressed (see [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF]), and the hydraulic pressure, tire-ground drag force and wheel bending load are introduced. By considering the interconnections between each element of the braking system, the equations of motions take the form:

Mẍ + C ẋ + Kx = F contact + F hydr + F drag + F load (1)
where ẍ, ẋ and x are respectively the 70-dimensional vectors of acceleration, velocity and displacement. F hydr is the vector force due to brake hydraulic pressure, F contact contains the linear and nonlinear contact force terms at the stator and rotor interface, F drag represents the tire-ground drag force, and F load the bending load of the wheel. M, C and K are respectively the structural mass, global damping and stiffness matrices.

Stability analysis

In a preliminary study, the stability analysis of the brake is performed. It is obtained for a small perturbation around the steady-state operating point by calculating the eigenvalues of the linearized system. The detailed methodology can be found in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF]. For λ a complex eigenvalue of the form λ = a + ıb, the stability criterion is: if a is negative or zero, the system is stable and no vibration occurs, and if a is positive then the system is unstable. The imaginary part b represents the angular pulsation (or the frequency by considering a factor of 2π) of the vibration. The plot of the angular pulsations versus the brake friction coefficient is given in Fig. 2. The graph is normalized towards the pulsation of the first instability ω 0,1 and towards the maximum friction coefficient µ max used in the study. The red plots define the appearances of instability. It appears that six of the fourteen firsts modes are involved in mode-coupling instabilities. The nature of each instability is determined by plotting its analytical shape. The mode shapes of all the instabilities can be found in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF]. The second instability will be the matter of the following study.

In [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] parametric studies are performed to evaluate the effects of damping or other parameters on the stability of the brake. The effects of damping on the second instability are summed up in Fig. 3, where the stability area is given when the damping ratio η 1 /η 2 between the two modes producing the instability varies. It turns out that there exists a value of the damping ratio for which the stability area is the widest. This optimum ratio is near 1, that means the optimum ratio corresponds to η 1 ≃ η 2 . In fact, in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] it is shown that the optimum damping ratio corresponds to a proportional damping configuration, given by η 1 /η 2 = ω 0,2 /ω 0,1 , where ω 0,1 and ω 0,2 are the pulsations at µ = 0 of the two modes producing the instability. But given that ω 0,1 and ω 0,2 are very close, the relation can be in first approximation expressed by η 1 ≈ η 2 , ie for a damping ratio close to 1. For values of the damping ratio far from the optimum value, the stability may be clearly affected, even leading to a more unstable system that the undamped system! That illustrates the phenomenon of "destabilization paradox" presented in introduction of the paper.

Introduction

Because damping appears to be a key factor, its effects on the dynamical behavior of the aircraft braking system must be examined. Determining the stable and unstable regions is only one aspect of the problem of friction-induced flutter instability. However, both the maximum amplitude of the instability and the time to reach the limit cycle are more significant design factors than the stability region.

The nonlinear equations of motion (1) can be integrated numerically to obtain time-history responses. For this study the case of the second instability (the second whirl) presented above is considered. For a value of µ/µ max between 0.045 and 0.115, this is the only instability that can occur (see Fig. 2), so it will simplify the analysis of the effects of parameters on the nonlinear behavior. In Fig. 4a the time-history response of the vertical deflection of the end of the axle is given. The simulation has been made with a friction coefficient µ/µ max of 0.075 and with an arbitrarily chosen damping configuration. One observes two different stages: first an initiation phase (I) where the oscillations increase around the static position to attain the stabilized phase (II) where the oscillations keep a constant amplitude. These two phases are well represented in Fig. 4b which shows the evolution in the coordinates system (displacement,velocity). The phase (I) is characterized by increasing spirals, and the phase (II) is the stationary limit cycle. With the application of a FFT transformation, it is verified that the signal contains only one vibration at a frequency very close of that given by the stability analysis with linear approximation (in a case of a nonlinear resonance, the frequency of the instability may be altered by nonlinear effects).

In order to understand the effects of structural damping on the oscillation amplitude and the time to reach the limit cycle, the next study is divided into two parts. In the first one, the effects of the friction coefficient and damping on the transient behavior are investigated. In the second one, the effects on the limit cycle and maximum amplitude are examined.

Investigations on the transient behavior

In this section the transient behavior corresponding to the initiation phase (I) is examined. The characterization of this phase is given by the response time: it represents the time to attain the limit cycle of the vibration. Here a criterion of 98% of the amplitude of the limit cycle is considered. In Fig. 5 is plotted the response time versus the brake friction coefficient. For a value of µ/µ max under 0.045, the system is stable and no vibration occurs. It can be seen that in this range of µ the response time is a decreasing function of the friction coefficient: the closer we are to the bifurcation point (where the two modes couple), the more slowly the vibration appears. But this conclusion requires more investigations. To further our understanding, the inverse of the real part of the mode becoming instable is displayed on the same figure. It is extracted from the eigenvalue analysis: if λ = a + ıb is the eigenvalue with a positive real part, then during the beginning of the phase (I) the displacement can be written in linear approximation by an exponential increase:

x = Ae λt = Ae at e ıbt ( 2
)
where A is a constant. Thus, the amplitude of the oscillation -of angular pulsation bis given by Ae at , with a time constant of 1/a. Let's note that the stability criterion appears here clearly: if a is positive the argument of the exponential is positive resulting in an increasing oscillation, and if a is negative or zero the amplitude of the oscillation is limited. In phase (II) the nonlinear effects induce the stabilization of the oscillation into the limit cycle. The Fig. 5 shows that a linear relation between the 98% response time and the inverse of the time constant of the exponential approximation can be established: t 98% ≈ 20 a . So the results of the stability analysis can be exploited in terms of time-history response, even if the numerical constant in the previous relation seems to be different for each instability. The stability analysis is very quick to obtain for a model of only 70 degrees of freedom. Then, in addition of the information obtained on the stability of the system, it also gives an idea of the temporal behavior. The second point is the evaluation of the effects of damping. Two approaches are undertaken. In the first case, proportional damping is considered: that means that the same damping coefficient η is introduced on the two coupling modes of the second whirl. In the second case, non-proportional damping is considered: the values of the damping coefficients η 1 and η 2 vary under the condition η 1 + η 2 = constant. Below are summarized the parameters of the two cases considered, both use a friction coefficient µ/µ max of 0.075. case #1 η 1 = η 2 , from 0.5% to 3% case #2 η 1 + η 2 = 4%, with η 1 /η 2 from 1/7 to 7

In the first case corresponding to proportional damping, the analysis is quite simple. Fig. 6a shows that the increase of proportional damping produces the slowdown of the vibration with a longer response time. Thus, from a damping of 0.5% to 3%, the time the vibration takes to attain the limit cycle goes from 0.8 seconde to 4.1 secondes, ie a ratio of more than 5. In the second case, the presence of an optimum is observed in Fig. 6b: a damping ratio close to 1 pushes back at best the vibration. That coincides with the results of the stability analysis stated in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] i.e. for this system proportional damping is the optimum configuration. In addition these two studies show once more the good relation between the response time and the real part of the unstable mode.

Investigations on the maximum amplitude of oscillations and limit cycle

This second section deals with the evaluation of the effects of friction coefficient and damping on the limit cycles and maximum amplitudes of oscillation in phase (II).

First, limit cycles are plotted for several friction coefficients. The chosen degrees of freedom are the horizontal and vertical deflections of the end of the axle, and the horizontal and vertical displacements of the middle of the brake rod. They are indeed elements in motion during the occurrence of the second (a) Proportional damping (b) Non-proportional damping Figure 6: Evolution of the 98% response time with variations of damping whirl. It may be observed that on average the increase of the friction coefficient leads to larger limit cycles. So, as expected, when moving away from the bifurcation point by increasing µ, the system becomes more unstable, which is illustrated by an increase of the oscillation amplitude and a faster to appear vibration.

Then, the effects of proportional and non-proportional damping are investigated thanks to the same two studies as those considered in the previous section. In Fig. 8 it is clear that the addition of damping is synonymous with the decrease of the size of the limit cycle. To have a more visual representation, in Fig. 10 is plotted the sum of the maximum amplitudes -of displacement or velocity-taken on the degrees of freedom of the brake rod on the one hand, and on those of the axle on the other hand. The evolution with the amount of damping shows a decrease of a factor of 11 between the configuration with 0.5% of damping and the configuration with 3% of damping. 3% can be assumed to be a realistic value for mechanical systems using damping elements, and an amount of 0.5% of damping can also be possible in the worst case. But the damping coefficient can hardly reach zero because of the damping properties of the materials and the always dissipative contacts between the elements of the structure. So the advantages of using damping in order to reduce friction-induced vibrations appear here clearly.

In the second study the effects of non-proportional damping are examined. Fig. 9 illustrates the evolution of the limit cycle with the damping ratio, and Fig. 11 the evolution of the maximum amplitudes of the brake rod and the axle. It can be observed that an optimum is obtained as well for the displacements as for the velocities. The best configuration (the one for which the amplitudes are minimal) is given by a damping ratio near 1, that confirms the results obtained in the section on the transient behavior. Considering a constant amount of damping, a change in the damping ratio can then result in oscillations almost two times bigger, which is absolutely not negligible. However, by placing damping elements only in the part of the system that vibrate the most, an asymmetry in the damping configuration can be obtained. This asymmetry could then be harmful because the mechanism could generate more vibration than the original undamped mechanism. A nonlinear model has been developed in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] for the stability analysis of mode-coupling instabilities induced by friction of an aircraft braking system. In this paper, the same model is used to evaluate the effects on the time-history response of the brake friction coefficient and the amount and ratio of damping. The first conclusion is that increasing the amount of damping leads to a vibration with a smaller amplitude and slower to develop. That points out the fact that damping can be very efficient to reduce or eliminate friction-induced vibrations. However, precautions have to be taken to avoid a misuse of damping. While in [START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF] it is shown that for this mechanism proportional damping is the best configuration for the stability of the system, the study presented leads to the same result about the dynamical behavior. At first sight it is not obvious for a nonlinear system because stability is obtained with linear approximation. Even if it may not be easy to design a system with an optimum damping distribution, in the light of the conclusions of this study, adding damping elements only in the part of the system that vibrate the most is an approach to avoid. In addition, the study shows that there exists a simple relation between the speed of increase of the vibration and the real part of the mode becoming unstable. Thus, the results of a stability analysis can be used in terms of temporal behavior, that is very interesting when time-history simulations are often time-consuming.

• • + • •η 1 /η 2 = 0.14, -+-η 1 /η 2 = 0.33, -+ -η 1 /η 2 = 0.6, • • • η 1 /η 2 = 1, -η 1 /η 2 = 1.67, --η 1 /η 2 = 3, -• -• η 1 /η 2 = 7) 0 0.5 1 1.5 2

Figure 1 :

 1 Figure 1: Aircraft braking system description

Figure 2 :

 2 Figure 2: Angular pulsations versus friction coefficient

Figure 3 :

 3 Figure 3: Stability (stable and unstable zones) versus µ/µ max and η 1 /η 2 for the second whirl instability

Figure 4 :

 4 Figure 4: Vertical displacement of the axle (µ/µ max = 0.075)

Figure 5 :

 5 Figure 5: Evolution of the 98% response time with variations of the friction coefficient

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Limit cycles -variations of µ(• • • µ/µ max = 0.05, --µ/µ max = 0.075, --µ/µ max = 0.1, • • + • •µ/µ max = 0.125, -+-µ/µ max = 0.15, -+ -µ/µ max = 0.175)

Figure 10 :

 10 Figure 10: Vertical displacement of the axle -proportional damping

Figure 11 :

 11 Figure 11: Vertical displacement of the axle -non-proportional damping