Denis Brizard
email: denis.brizard@ec-lyon.fr

Olivier Chiello
email: olivier.chiello@inrets.fr

Jean-Jacques Sinou

Xavier Lorang

J.-J Sinou

Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact

come

Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact

n, t, b normal, tangential and radial directions at a particular contact point P n , P t projection matrices on the relative displacements between disc and pads along the normal and tangential directions r n , r t normal and tangential reaction forces at the contact nodes T c basis for the kernel of the normal contact constraints P n T r some reduced basis or other M, K mass and stiffness matrix of the uncoupled system Mc , Kc mass and stiffness matrices of the Coupled system with friction M µ , K µ mass and stiffness matrices taking into account the effects of the friction forces C b radial dissipation effect damping matrix K c stiffness matrix of the frictionless coupled system µ friction coefficient ω i , R i , ξ i imaginary part (frequency), real part and divergence rate associated with mode i ωi , ξi circular frequency and divergence rate of mode i of reference solution Φ, Ψ matrix of normal modes, matrix of static modes n m , n c number of retained normal Modes, number of Contact nodes Superscripts *, T conjugate and transpose of a matrix or a vector ˆrelated to reference solution D, P related to Disc and Pins Subscripts f , s related to free boundary normal modes and sliding boundary normal modes an, cn related to Normal Attachment modes and Normal Constrained modes at, aµ related to Tangential Attachment modes and rubbing Attachment modes R related to the normal modes of the frictionless coupled structure cat related to the Tangential Attachment modes of the Coupled structure

Introduction

The noise generated by vehicles brakes is one of the most difficult problems with brake systems. Research for predicting and removing brake noise has been regularly performed for many years, especially for the automotive industry [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal: a review[END_REF][START_REF] Crolla | Brake noise and vibration -State of the art[END_REF][START_REF] Yang | Brake vibration and noise: reviews, comments and proposals[END_REF][START_REF] Papinniemi | Brake squeal: a literature review[END_REF][START_REF] Kinkaid | Automotive disc brake squeal[END_REF]. Among the different encountered brake noises, the squeal noise is characterized by very annoying spectra with one or several pure frequencies in the medium or high frequency range (1-15 kHz). Despite great progress in the modeling of brake squeal, many of the noise abatement approaches are still individual and empirical. It is therefore an important issue to be able to predict this noise.

Friction between two structures in contact can generate many kinds of vibrations and noise [START_REF] Akay | Acoustics of friction[END_REF]. Some of them are forced vibrations, for instance due to the macroscopic roughness of the surfaces in contact. Others are parametric vibrations due to the moving loads. However, a common view is that squealing is, above all, the result of selfsustained vibrations of the brake system [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal, and chaos part II: Dynamics and modelling[END_REF]. There are mainly two types of instability that can lead to self-sustained vibrations. Both are due to the non-conservative nature of the friction forces. The first highlighted cause of instability is the possible decrease of the friction coefficient with the sliding speed [START_REF] Hinrichs | On the modelling of friction oscillators[END_REF][START_REF] Bengisu | Stability of friction-induced vibrations in multi-degree-of-freedom systems[END_REF]. However, most brake squeal models are based on a second instability mechanism which is often called 'mode coupling' [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF][START_REF] Hervé | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF][START_REF] Ouyang | Unstable travelling waves in the friction-induced vibration of discs[END_REF][START_REF] Lorang | TGV disc brake squeal[END_REF]. Geometric instabilities or kinematically constrained instabilities are the other employed names for this mechanism. This may occur with a constant friction coefficient, as defined by Coulomb's law, which couples normal and tangential forces at the interface.

In order to determine the stability of a dynamical system, a classical method consists in solving the eigenvalue problem corresponding to the system linearized around the equilibrium [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF]. This simple technique provides complex modes and complex eigenvalues. Eigenvalues with positive real parts are characteristics of unstable equilibrium points which diverge from steady sliding and lead to self-sustained vibrations. Many studies are limited to this approach and do not explicitly calculate the self-sustained vibrations. Indeed, a pure non-linear transient or stationary analysis is necessary to compute the limit cycles the problem tends to [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].

When finite element models are considered, computing the stability requires to solve a non-symmetrical eigenvalue problem with a great number of degrees of freedom. Such a problem can be solved by using widespread iterative [START_REF] Sorensen | Numerical methods for large eigenvalue problems[END_REF] or less known techniques such as the residue iterative method [START_REF] Bobillot | Méthode de réduction pour le recalage[END_REF]. In any case, the use of iterative methods is highly expensive. To overcome this computational problem, a common technique is to expand the dynamic equations on reduced real bases. Practically, bases composed of frictionless (real) coupled modes or merely free components modes are often used. However, the question of errors provided by the reduction is rarely investigated.

In this context, Component Mode Synthesis (CMS) [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF][START_REF] Mac | A hybrid method of component mode synthesis[END_REF][START_REF] Rubin | Improved component mode synthesis for structural dynamic analysis[END_REF][START_REF] Craig | On the use of attachment modes in substructure coupling for dynamic analysis[END_REF][START_REF] Craig | A review of time-domain and frequency-domain component-mode synthesis[END_REF] appears to be a logical issue for the dynamical study of a brake mechanism composed of many substructures. It offers the advantage of computing component modes separately, thus requiring less computational efforts than directly calculating the global modes. A key point in CMS techniques is that the reduced bases not only contain components' normal modes but also some static boundary modes, such as attachments modes in free-interface methods or constrained modes in fixed-interface ones [START_REF] Craig | A review of time-domain and frequency-domain component-mode synthesis[END_REF]. These modes allow statically complete bases to be built which greatly improve the results.

In this paper, the performances of reduced bases are investigated in the case of a simple disc/pads system. Bases composed of real coupled modes and bases provided by CMS techniques are tested. An enrichment of these bases is proposed in order to improve the precision of the results. The idea behind this enrichment is to build bases statically complete even when friction forces are taken into account. In particular, new rubbing attachments modes are proposed to adapt free-interface CMS techniques to frictional contact. In addition, when real coupled modes are used, it is suggested to complete the basis by the static response to a first-order approximation of the friction forces.

In section 2, the formulation of the problem is described, based on the previous work initially proposed by Moirot [START_REF] Moirot | Etude de la stabilité d'un équilibre de Coulomb -Application au crissement des freins à disque (Stability study of a Coulomb equilibrium -Implementation on disc brake squeal)[END_REF] and used for instance by Lorang [START_REF] Lorang | TGV disc brake squeal[END_REF]. Reduction principle and errors criteria which are employed to evaluate the reduction approximation are given in section 3. Section 4 introduces the disc brake model. Components modes and real coupled modes are analyzed and a preliminary stability study is performed. Classical CMS bases are first applied on the frictionless problem in section 5. Finally, in section 6, enrichments are introduced for the problem with friction and performances of classical and enriched bases are compared and discussed.

Formulation of the problem

In this section, the discretized equations governing the stability of the steady sliding of a disc/pads system are derived. The main assumptions are the linear elasticity in the volumes of the structures combined with unilateral contact and Coulomb friction at the interface. At each point of the interface, a local cartesian coordinate system (n, t, b) is defined where n is the direction normal to the disc, t is the direction of the steady sliding and b is the radial direction defined as b = n ∧ t (c.f. figure 1).

Continuous non-linear form

The rotation speed of the disc is supposed to be constant and sufficiently low so that the gyroscopic terms and the volume forces induced by rotation may be neglected. By using an eulerian description, the displacements field u and the stresses tensor σ in the volumes D and P of the disc and the pads verify:

σ = A : ∇u, ∇ • σ = ρü ∀x ∈ D, P (1)
where A is Hook's tensor.

On the surface S c potentially in contact, the disc and the pads are coupled with conditions of unilateral contact with Coulomb friction, which may be expressed by:

R n ≤ 0, g ≤ 0, gR n = 0 ∀x ∈ S c R t ≤ -µR n , w 0 ⇒ R t = µR n w w ∀x ∈ S c (2)
In these equations, g = g 0 + u n denotes the normal gap between the disc and the pads where g 0 is the initial gap and u n is the relative normal dynamic displacement. In addition, w = w 0 t + ut denotes the sliding velocity between the disc and the pads where w 0 is the stationary sliding velocity due to the disc rotation and ut is the relative tangential dynamic velocity. Finally, R n and R t represent the normal and tangential contact forces between the disc and the pads.

The friction coefficient µ is supposed to be constant.

Continuous linear form

If the rotation speed is not equal to zero, the previous equations may have a quasi-static solution, or steady sliding, for which u = 0 (the existence and the uniqueness of this equilibrium is not discussed here). If initial conditions differ from the steady sliding, two cases may be distinguished: the stable -non-squealing -case if the solution converges to the equilibrium, and the unstable -squealing -case if the solution grows, leading to self-sustained vibrations. The analysis of the stability is generally achieved by a linearization method. Considering small regular sliding perturbations that do not break the contact, the contact forces may be linearized. On the contact surface S 0 in the quasi-static solution, the perturbation verifies the following equations:

u n = 0 ∀x ∈ S 0 R t = µR n t -c b ub b ∀x ∈ S 0 (3)
In these equations, ub = ut .b denotes the relative tangential velocity between the disc and the pads in the direction b and c b = R n0 /w 0 where R n0 is the quasi-static normal contact force. In the expression of the linearized friction force, the first term (along t) is classical and comes from the linearization of the module of the friction force. The second term (along b) is less well-known and is due to the linearization of its direction ; it only appears in three dimensional models and acts as a viscous contact damping (it is also referred to as radial dissipation effect by some authors [START_REF] Hochlenert | Friction induced vibrations in moving continua and their application to brake squeal[END_REF]).

Discrete linear form

By using a finite element method with conforming linear meshes at the interface, the previous linear equations ( 1) and (3) may be written in a discrete form as follows :

M ü + C b u + Ku = P T n r n + P T t r t P n u = 0 r t = µr n (4)
where M and K denote the mass and stiffness matrices of the uncoupled system, C b is the damping matrix associated to the viscous contact damping c b and u represents the corresponding displacements vector. In addition, r n and r t denote the normal and tangential reaction forces at the contact nodes along the directions n and t whereas P n and P t are the projection matrices on the relative displacements between the disc and the pads along the same directions. After some algebraic manipulations, all the linear contact constraints may be eliminated and system (4) may be written in a more compact form :

u = T c u c Mc üc + Cc uc + Kc u c = 0 (5)
where

Mc , = T T c MT c Kc = T T c KT c Cc = T T c C b T c (6)
and

M = M + M µ K = K + K µ (7)
In the above expressions, T c is a basis for the kernel of the normal contact constraints P n , M µ and K µ are nonsymmetrical matrices taking into account the effects of the friction forces and given by:

M µ = µ 2 P T t P n M K µ = µ 2 P T t P n K (8)

Complex eigenvalue analysis

The stability of the steady sliding may then be deduced from a complex eigenvalue analysis of the system:

λ 2 Mc + λ Cc + Kc φ c = 0 (9)
providing a set of complex modes φ i = T c φ ci and complex eigenvalues λ i = R i + iω i in the frequency range of interest. The perturbation may be approximated by a truncated superposition of these modes, oscillating at their natural pulsation ω i but with divergent or vanishing envelopes controlled by R i and initial complex amplitudes a i = |a i | exp(ϕ i ) provided by the initial conditions :

u Re        i a i φ i exp(λ i t)        i |a i | exp (R i t) (Re (φ i ) cos (ω i t + ϕ i ) -Im (φ i ) sin (ω i t + ϕ i )) (10)
If all the modes are stable (i.e. ∀i R i < 0) the perturbation vanishes. On the other hand, if a complex mode is unstable (i.e. ∃i R i > 0), the perturbation grows, leading to self-sustained vibrations. This case may happen with friction forces because of the presence of the non-symmetric matrices M µ and K µ . For each mode, a divergence rate may also be defined, physically equivalent to a negative damping ratio :

ξ i = R i ω i (11)
R i (unit s -1 ) corresponds to the geometric divergence of the envelope of mode i, whereas divergence rate ξ i (unit rad -1 ) corresponds to that geometric divergence in phase unit.

Reduction of the problem

Principle

Solving a large non-symmetric eigenvalue problem like system (9) is very expensive from a computational point of view. Iterative algorithms are available but computation times are often prohibitive to conducting parameter studies.

Consequently, Ritz reduction methods are often used to solve the problem. This consists in searching for the best approximated solution in a reduced subspace characterized by a basis T r , by solving the reduced eigenvalue problem:

λ 2 Mr + λ Cr + Kr φ r = 0 ( 12
)
where

Mr = T T r MT r Kr = T T r KT r Cr = T T r C b T r (13)
providing a set of approximated complex modes { φi } = T r φ ri and approximated complex eigenvalues λi = Ri + i ωi in the frequency range of interest. Obviously, to be compatible, vectors in the reduced subspace must at least verify the normal contact constraints, i.e. P n T r = 0.

A widely used method is to choose the real coupled modes (without friction) of the structure as a reduced basis. Indeed, these modes are calculated by solving a real symmetric eigenvalue problem, for which very efficient iterative algorithms are available. Another used method is to build a reduced basis from component modes, according to the so-called Component Mode Synthesis. These bases may give good results in some cases. However, as shown in the last section of the paper, a better approximation is obtained by using bases which take into account particular static responses to friction forces.

Quality of the approximation

First, the criteria employed to evaluate the quality of the approximation due to the reduction are presented below.

Error on the residues

A first criterion is defined by considering the force residue of the non-reduced eigenvalue equation for a couple of approximated mode φi and approximated eigenvalue λi :

R F i = T T c λ2 i M + λi C b + K φi ( 14
)
where the matrix T c is needed here to eliminate the normal reaction forces at the interface. The method proposed by Bobillot is then used [START_REF] Bobillot | Méthode de réduction pour le recalage[END_REF]: in order to obtain an energy indicator, the force residue is first converted into a displacement residue by using the stiffness matrix of the frictionless coupled system K c = T T c KT c :

R U i = K -1 c R F i (15)
Finally, the residue error is defined as the ratio between the deformation energy of the residue and the deformation energy of the considered mode :

ε res i = R U i * T R F i φi * T K φi (16)

Error on the mode shapes

Another criterion is defined by comparing the approximated complex modes with the exact complex modes. For a given exact mode φ i , the first step is to identify the corresponding approximated complex mode. This is achieved by reordering the approximated base such that: where MAC-M is the complex mass-weighted Modal Assurance Criterion defined by:

MAC-M φ i , φi = max j MAC-M φ i , φ j (17)
MAC-M φ i , φ j = φ * T i M φ j 2 φ * T i Mφ i φ j * T M φ j (18)
This criterion was already defined for real modes as the NOC [START_REF] Morales | Comments on the mac and the nco, and a linear modal correlation coefficient[END_REF]. Thanks to this reordering, φi is the closest approximated mode to φ i by considering the classical complex norm and scalar product associated with the frictionless mass matrix. The error on the mode shapes is then defined as:

ε mac i = 1 -MAC-M φ i , φi (19)
No criterion is defined on the extra diagonal terms of the MAC-M since complex modes are not supposed to be orthogonal with each other. Indeed, no hermitian scalar product may be associated with the complex eigenvalue problem because of the non-symmetry induced by the friction forces.

Error on the frequencies and divergence rates

The matching of the approximated and exact bases now allows the calculation of some errors on frequencies:

ε frq i = |ω i -ωi | ω i (20)
and also on divergence rates:

ε div i = |ξ i -ξi | (21)
where ξ i is the exact divergence rate and ξi the approximated divergence rate. This last criterion is defined as an absolute error because the divergence rates are relative indicators themselves. Practically, in order to make any conclusions on stability, this error has to be compared with the precision of a damping ratio measurement.

The disc/pads system

In this section, the model used to evaluate the performance of some reduced bases is described. This simple model (cf. figure 2) is derived from a more complex one which was developed to find instabilities in a TGV (French highspeed train) braking system [START_REF] Lorang | TGV disc brake squeal[END_REF][START_REF] Chiello | Experimental characterization and numerical modelling of the squeal noise generated by railway disc brakes[END_REF][START_REF] Chiello | Numerical investigations into the squeal propensity of a railway disc brake[END_REF][START_REF] Lorang | Stability and transient analysis in the modelling of railway disc brake squeal[END_REF]. The main characteristic of this system is the geometry of the pads and in particular the friction material which is distributed into a set of small cylindrical pins. Compared with the original model, several simplifications have been performed in order to reduce its size. Firstly, the flexibility of the pads' supporting structures (backplate and support) is neglected: a fixed boundary condition is assumed for the external face of the cylindrical pins. Secondly, only 6 pins applied on either sides of the disc are taken into account (vs 18 for the original model). Lastly, a rougher mesh is used since the frequency range of interest is restrained here to 0 -5000 Hz which leads to a total number of degrees of freedom n = 30942 and a number of contact nodes n c = 228. This

Component modes

In a first step, the modes of each component are determined separately in order to go further into their own dynamics. For the disc with free boundary conditions on the contact surface, three mode types are distinguished:

• the axial modes Am-n for which axial deformations are dominant (out-of-plane, flexion);

• the radial modes Rm-n for which radial deformations are dominant (in-plane);

• the circumferential modes Cm-n for which circumferential deformations are dominant (in-plane).

In addition, each mode is characterized by a number of nodal circles m and a number of nodal diameters n. Because of the axi-symmetry of the disc, all modes with one nodal diameter or more are double, which means that there are two modes at the same frequency with the same shape but shifted by an angle equal to π/(2n). Some examples of the different encountered modes are represented on figure 3.

Pins modes are computed for 'sliding' boundary conditions on the contact surface that are fixed in the normal direction and free in the tangential one. With these conditions, the dynamic behavior of the pins is rather close to the one which takes place in coupled conditions since the disc is much more rigid than the pins. Only three modes per pin are found below 5 kHz: a double flexural mode F1 and a torsional mode T1 (cf. figure 4).

Real coupled modes

The second step is the determination of the real coupled modes (without friction i.e. M µ = K µ = C b = 0) for which the disc and the pins are coupled only through normal degrees of freedom on the contact surface. They are computed here by using a classical IRA/Sorensen iterative algorithm [START_REF] Lehoucq | ARPACK Users' Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods[END_REF]. In the frequency range above 1500 Hz, these coupled modes are easy to describe. One may distinguish the modes controlled by the disc from the modes controlled by the The first ones are close to the modes of the free disc presented above but with pins deformations following the disc motion, whereas the other ones are very similar to the sliding pin modes F1 and T1 but with a little disc contribution (cf. figure 5). In the low-frequency range, the higher stiffness brought by the pins tends to modify more strongly the free disc modes and some of them are harder to describe. However, a designation based on components modes is kept for most of the coupled modes in all the frequency range (cf. table 1).

                  
It is important to notice that, as the coupling breaks the symmetry of the components, disc double modes are split and pins multiples modes are dispersed. In particular, the difference between the frequencies of two disc modes of the same type (calledand +) may be explained by the difference of global stiffness brought by the pins according to the position of the disc nodal diameters. This is very pronounced for the axial disc modes and flexural pins modes. On the other hand, coupling has very little effect or no effect at all on some component modes with small normal contact amplitude such as the in-plane disc modes (circumferential and radial) and torsional pins modes.

Complex modes

The stability of the system is studied here for 90 values of the friction coefficient µ, varying from 0 to 0.45, by using the residue iterative method [START_REF] Bobillot | Méthode de réduction pour le recalage[END_REF]. Two cases are considered by taking the damping matrix C b into account or not ; note that no structural damping is taken into account in this paper. This matrix is the only one which depends on the equilibrium and thus is sensitive to the braking conditions. In this study, it is determined for a rotation speed equal to 4 rad/s and a static global load of 1000 N on each side of the disc. Neglecting it amounts to considering a high speed/low effort configuration. In both cases, the obtained complex modes are close to the real coupled modes in terms of frequencies and shapes except that some modes become unstable whereas others are damped.

Stability for µ = 0.4

When matrix C b is neglected, the system has four unstable modes, two disc modes and two pins modes (cf. table 1). The main characteristic of the unstable modes is that they are not stationary: they rather look like rotating waves. Moreover, they work in pair with a similar stable mode, rotating the opposite way with a divergence rate of the same magnitude but the opposite sign. Adding damping matrix C b , while there is no structural damping, has a 6) by reducing the divergence rates. As a consequence, the second slightly unstable pins mode becomes stable and many pins modes turn to highly negative divergence rates. It has been demonstrated that increasing damping does not always have a stabilizing effect [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF]. The system is stabilized by a "lowering effect" when the two modes involved in the coalescence are equally damped. On the contrary, both "lowering" and "smoothing" effects take place when the two modes are not equally damped ; in the case of a prevailing "smoothing effect" the system may be unintuitively destabilized by added damping. However the influence of damping upon stability is beyond the scope of this paper.

Bifurcation

Observing the evolution of the unstable modes with varying friction coefficient is useful to clarify the origin of instabilities. When neglecting C b , the coalescence of the disc modes is easy to observe (cf. figure 7): the frequencies of the two similar disc modes, distinct at µ = 0, become equal with increasing friction coefficient. There are then two modes at the same frequency, one being stable and the other unstable. The same coalescence phenomenon occurs for pins modes but the higher modal density makes it a bit more difficult to see (cf. figure 8). With increasing friction coefficient, an initially unstable pin mode restabilizes for µ = 1.8. When adding C b , rather similar results can be observed (cf. figure 7 and figure 9). Computing these graphs is extremely expensive: a complex eigenvalue problem has to be iteratively solved for each friction coefficient step. A reduced basis could therefore be useful to reduce the computation times.

Reduced bases for the frictionless problem

Before using reduced bases for the stability study, component mode synthesis methods (CMS) are tested on the problem without friction. The objective is to evaluate the performance of some bases constructed from the component modes for the calculation of the real coupled modes. This is a classical problem but it is presented here as a preliminary study and also as an introduction to the formalism used further in the stability study of section 6.

Description of the bases

Reduction of the components

For the reduction of the components, two classical CMS bases are used. The first one is a combination of dynamic free-interface normal modes and static boundary attachments modes [START_REF] Mac | A hybrid method of component mode synthesis[END_REF][START_REF] Rubin | Improved component mode synthesis for structural dynamic analysis[END_REF][START_REF] Craig | On the use of attachment modes in substructure coupling for dynamic analysis[END_REF] whereas the second one uses fixedinterface dynamic normal modes plus static boundary constrained modes (original Craig-Bampton method [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF]). In both cases, static boundary modes are added to dynamic normal modes in order to make the bases statically complete for the coupled system. In our problem, two components are considered (the disc and the whole set of pins) and the interface is restrained to normal degrees of freedom of the contact surface. In addition, since the disc and the pins have no rigid body modes, the definitions of the static boundary modes do not present any particular problem. If one then defines:

• the matrix Φ D f (resp. Φ P f ) of the n D m (resp. n P m ) first normal modes of the disc (resp. pins) with free boundary conditions on contact nodes,

• the matrix Φ D s (resp. Φ P s ) of the n D m (resp. n P m ) first normal modes of the disc (resp. pins) with sliding boundary conditions on contact nodes i.e. fixed in the normal direction and free in the tangential one,

• the matrix Ψ D an (resp. Ψ P an ) of the n c normal attachments modes of the disc (resp. pins) i.e. the static responses of the disc (resp. pins) to imposed unit normal forces on contact nodes,

• the matrix Ψ D cn (resp. Ψ P cn ) of the n c normal constrained modes of the disc (resp. pins) i.e. the static responses of the disc (resp. pins) to imposed unit normal displacements on contact nodes, the disc (resp. pins) free-interface complete base may be expressed as Φ D f Ψ D an (resp. Φ P f Ψ P an ) and the disc (resp. pins) fixed-interface complete base may be written Φ D s Ψ D cn (resp. Φ P s Ψ P cn ). In our notations, all the modes are defined for the same set of degrees of freedom corresponding to the global displacement vector u.

Connection of the components

Considering the two kinds of component bases presented above, four bases may be defined for the whole structure according to the choice made for each component. In this study, only three of them are tested (see also table 2):

• the "free" one combines free-interface bases for both the disc and the pins ; in this particular case, a global set of static boundary attachments modes (Ψ D an -Ψ P an ) is used rather than one for each component ; this is a frequently used free-interface CMS technique : it amounts to adding constraints on the generalized coordinates (representing static boundary forces) associated to the attachments modes [START_REF] Craig | On the use of attachment modes in substructure coupling for dynamic analysis[END_REF][START_REF] Craig | A review of time-domain and frequency-domain component-mode synthesis[END_REF];

• the "fixed" one combines fixed-interface bases for both the disc and the pins; • the "mixed" one combines a free-interface basis for the disc and a fixed-interface basis for the pins.

Method Reduced basis T r0 size of T r0 size of T r free Φ D f Φ P f (Ψ D an -Ψ P an ) n D m + n P m + n c n D m + n P m fixed Φ D s Φ P s Ψ D cn Ψ P cn n D m + n P m + 2n c n D m + n P m + n c mixed Φ D f Φ P s Ψ D an Ψ P cn n D m + n P m + 2n c n D m + n P m + n c
Actually, in the mixed basis, the component normal modes presented in section 4 are used. It is thus expected to exhibit good results since the real coupled modes are close to these component modes. On the contrary, a mixed basis with a fixed-interface basis for the disc and a free-interface basis for the pins would be of poor interest since, in this case, component modes would be very far from the real coupled modes. Before using the basis to solve the eigenvalue problem, the connection of the components is achieved by forcing the base to satisfy the normal contact constraints P n T r = 0. For an initial unconstrained reduced base T r0 , this is done by taking a final constrained reduced base T r such that:

T r = T r0 T cr ( 22
)
where T cr is a basis for the kernel of P n T r0 . It generally leads to the elimination of n c vectors in the basis which results in the bases sizes given in table 2.

Results

In order to obtain good results with CMS, it is well known that component normal modes need to be computed in an extended frequency range compared with the one of interest. The three methods have thus been tested for three sets of normal modes corresponding to a maximum frequency equal to 1.2, 1.5 and 2 times the limit frequency of the study i.e. 6, 7.5 and 10 kHz. The numbers of disc and pins modes corresponding to these frequencies are given in table 3. The errors due to the reduction are given on figures 10, 11 and 12. The error criteria on residues ε res i , frequencies ε frq i and mode shapes ε mac i defined previously are used (which corresponds respectively to the res, frq and mac labels on the x axis of each figure). For the errors on frequencies and shapes, the reference modes are the real coupled modes presented in section 4 and computed by using a classical IRA/Sorensen algorithm. Only the maximum errors for each kind of mode (disc modes, F1 pin modes and T1 pin modes) are given.

For all tested bases, errors on torsional T1 pin modes are extremely small (out of scale on the figures) except for criterion ε mac i . Actually, these modes have the same frequency and are not coupled with the disc: these are purely multiple modes. As a consequence, all linear combinations of them are solutions of the eigenvalue problem and a criterion based on individual shapes comparison is not pertinent in this particular case. For the other kinds of modes, results are coherent: errors corresponding to the three criteria move together in the same direction with the same order of magnitude. However, errors ε mac i tend to reduce with more difficulty than the others. The free basis gives rather good results, even in the 6 kHz case, especially for the disc modes. This is not surprising since this is the disc which controls the dynamics of the whole system. When component modes up to 7.5 kHz or 10 kHz are kept, the errors are even divided by 100 for both the disc and pins modes and the free method reaches Error criterion free 10kHz fixed 10kHz mixed 10kHz

Figure 13: Global comparison of the performance of the CMS bases for the frictionless problem its optimum. Conversely, the fixed basis is not adapted for a good description of the disc modes but gives acceptable results for the pins modes. However, keeping more components modes does not improve the performances of the basis in a significant way. As expected, the mixed basis gives very small errors for both the disc and the pins with a reduced number of component modes. However, as for the fixed method, the errors are not significantly reduced when the number of component modes increases. A global comparison of the three bases is given on figure 13. For this comparison, the maximum error on both the disc and F1 pins modes is taken into account and components modes up to 10 kHz are kept. It is now clear that, despite the initial good performances of the mixed basis, the free one gives the best results when it has converged: the dynamics of the coupled system is better described when a sufficient number of free pins modes are taken into account.

Reduced bases for the stability study

In this section, the use of reduced bases for the stability analysis is investigated. The friction is taken into account in the eigenvalue problem leading to complex eigenvalues and complex modes. Two kinds of bases are tested and compared together:

• bases built from the components modes (CMS),

• bases built from the real coupled modes.

Enrichment of CMS bases

The CMS bases defined in the previous section and adapted to the frictionless problem may be used for the stability study. However, since friction is added, some forces apply on the tangential contact degrees of freedom in the coupled system. Therefore, it may be interesting to add some static responses to these forces in the bases through the use of tangential attachment modes. To that end, we define the matrix Ψ D at (resp. Ψ P at ) of the n c tangential attachments modes of the disc (resp. pins) ; it corresponds to the static responses of the disc (resp. pins) to imposed unit tangential forces on contact nodes in the direction t. Thus, the enriched CMS bases given in table 4 may be constructed. An even more adapted approach is to include in the basis the static responses to normal and friction forces applied together on the components. This is done by defining some new rubbing attachment modes for the disc and the pins as:

Ψ D aµ = Ψ D an + µΨ D at Ψ P aµ = Ψ P an + µΨ P at (23)

Method

Reduced basis T r0 size of T r0 size of T r Free + t.a.m.

Φ D f Φ P f (Ψ D an -Ψ P an )(Ψ D at -Ψ P at ) n D m + n P m + 2n c n D m + n P m + n c Free with r.a.m. Φ D f Φ P f (Ψ D aµ -Ψ P aµ ) n D m + n P m + n c n D m + n P m Fixed + t.a.m. Φ D s Φ P s Ψ D cn Ψ P cn (Ψ D at -Ψ P at ) n D m + n P m + 3n c n D m + n P m + 2n c Mixed + t.a.m. Φ D f Φ P s Ψ D an Ψ P cn (Ψ D at -Ψ P at ) n D m + n P m + 3n c n D m + n P m + 2n c
Method Reduced basis T r size of T r Real Φ R n R Real + t.a.m. Φ R (Ψ D cat -Ψ P cat ) n R + n c Real + 1 st -ord. Φ R Ψ Rµ 2n R
Table 5: Bases constructed from the real coupled modes for the stability study (t.a.m.: tangential attachment modes, 1 st -ord.: 1 st -order approximation of the response to frictional forces)

However, this enrichment technique may only apply when normal attachments modes are already present in the initial basis, which is the case for the free and the mixed CMS bases. In this paper, it has only been tested for the enrichment of the free CMS basis as it is summarized in table 4.

Enrichment of the real coupled modes

As previously mentioned, a widely used method is to choose the real coupled modes of the structure as a reduced basis. It is pertinent since efficient iterative algorithms are available for the resolution of the symmetric real problem. Thus, we define the matrix Φ R of the first n R normal modes of the frictionless coupled structure. Although this basis may give rather good results, it would be interesting to complete it by adding some static response to the friction forces, similar to the enrichment of the CMS bases. This is done by using the matrix Ψ D cat (resp. Ψ P cat ) of the n c tangential attachments modes of the coupled structure (with contact constraints but no friction). It corresponds to the static responses of the coupled structure to imposed unit tangential forces on disc (resp. pins) contact nodes in the direction t. With this technique, the basis is statically complete since the responses to every possible friction forces are included.

Another possibility is to carry out a first-order approximation of the friction forces and the corresponding static response. Although real coupled modes are obtained without friction, it is possible to calculate the static part of the friction forces which would be induced by the displacement field due to these modes:

F (Φ R ) = K µ Φ R (24)
With the use of the stiffness matrix K c of the frictionless coupled system, the corresponding static response is obtained:

Ψ Rµ (Φ R ) = T c K -1 c T T c K µ Φ R (25)
Matrix Ψ Rµ may then be used to complete the real coupled modes. This enrichment method finds its origin in the work of Plouin [START_REF] Plouin | Steel/viscoelastic/steel sandwich shells, computational methods and experimental validations[END_REF] where it has been proposed in a similar form for viscoelastic problems. The two resulting enriched bases and their size are summarized in table 5. The errors due to the reduction are given on figures 14, 15 and 16 for the CMS bases and on figure 17 for the bases built from real coupled modes. The error criteria on residues ε res i , frequencies ε frq i , divergence rates ε div i , and shapes ε mac i previously defined are used (which corresponds respectively to the res, frq , div and mac labels on the x axis of each figure). For the errors on frequencies, divergence rates and shapes, the reference modes are the complex modes presented in section 4 and computed by using a residue iterative algorithm. In the CMS bases, the maximal number of component modes is taken into account (see table 3 at 10 kHz). In the bases constructed from the real coupled modes, a number of normal modes n R = 66 is generally used which corresponds to a truncation frequency of 6 kHz. A test with a more extended real base with n R = 130 modes, i.e. a truncation frequency of 10 kHz, is also performed for comparison.

As for the frictionless case, the maximum errors for each kind of mode are given but unstable and stable pins modes are distinguished instead of F1 and T1 pin modes. Indeed, the problem of pertinence of the shape criterion ε mac i for multiple modes is also encountered here for some sets of stable F1 pin modes which are very close to each other. They are grouped together with the T1 pins modes which are still not affected by the contact, even with friction (cf. table 1). Except this particular problem for some quasi-multiple stable pins modes, the different errors criteria give rather coherent results. However, errors on mode shapes tend to decrease with even more difficulty than in the frictionless case and errors on divergence rates seem to be very sensitive in some cases.

The results obtained when using the CMS or real classical bases are acceptable, but it is clear that adding new static modes adapted to friction improves their performances. This is particularly true for the bases built from free components modes or real coupled modes (cf. figure 14 and17) which share a very good initial description of the dynamics of the frictionless system. With the newly introduced static modes, all errors are reduced by a factor comprised between 10 and 100 both for the disc modes and the unstable pins modes. On the other hand, the extension As expected, the fixed CMS basis, which was not efficient for the frictionless problem, still exhibits poor performances for the stability study even if enriched (cf. figure 15). The mixed basis gives intermediate results (cf. figure 16). Without tangential attachments modes, it provides surprisingly low errors for the divergence rates of disc modes. As a consequence, when it is enriched, the errors are reduced except for these divergence rates. However, when considering the maximal error for both disc and pins modes, the enriched basis is at least as effective as the initial one whatever the criterion.

A global comparison of the enriched bases is given on figure 18. For this comparison, the maximum error on both the disc and unstable pins modes is taken into account. It appears that the performances of the bases built from free components modes and real coupled modes are similar. From a practical point of view, these are the only bases which provide a maximal error on divergence rates lower than 0.1% (10 -3 in absolute value on figure 18) which seems to be necessary to correctly determine the stability of a lightly damped structure. In table 6, the size n r of these four bases, as well as the total computational time t r (including the building of the base) needed to calculate the complex modes are given. Relative indicators are also calculated:

• the size gain of the basis G n = n n r where n = 30942 is the number of degrees of freedom of the whole structure,

• the computational time gain of the basis G t = t t r where t = 3242 s is the time needed to solve the complex eigenvalue problem without reduction.

For these four bases, the order of magnitude of the size gain is 100 whereas the cpu time gains are all greater than 20, which are both considerable values. The use of real coupled modes combined with a 1 st -order approximation of the response to frictional forces is a very interesting option, providing the highest gains both in size and cpu time. If, for some external constraints like experimental updating or damping modeling, components modes are preferred, the use of rubbing attachments modes seems to be a good technique to complete the free-interface components normal modes. Indeed, in this case, the size of the basis depends only on the number of components normal modes and not on the number of contact nodes (cf. table 4). Conversely, the two methods using attachments modes provide bases with higher sizes. However, their time gains are rather similar and for systems where the number of contact nodes is small compared with the number of modes in the frequency range of interest, the use of tangential attachment modes may be a better option to complete either the components modes or the real coupled modes.

Method

Conclusion

The complex eigenvalue analysis is a widely used technique to investigate the stability of a dynamical system with frictional contact. In the case of brake systems, it is the most frequently employed method to study the propensity of the brake to generate squeal noise. When iterative solvers are used, the complex modes and eigenvalues are calculated with good precision. In practice, reduced real bases are often used in order to save computational times. In this case, great attention should be focused on the choice of the basis. Indeed, it has been shown that classical CMS techniques or bases composed of real coupled modes may lead to significant errors in terms of frequency, divergence rates and shapes. In particular, errors on divergence rates greater than 0.1% lead to an erroneous prediction of the stability of a lightly damped system which may have important practical consequences.

In this work, some enriched bases are proposed in order to improve the precision on the calculated complex modes and eigenvalues. When bases built from component modes are used, it is suggested to add tangential attachment modes to classical CMS bases or to use the new rubbing attachment modes defined in the paper. In the case of the disc/pads system studied here, such enrichments on free-interface CMS bases allow a significant decrease of the errors on frequencies, divergences rates and shapes by a factor comprised between 10 and 100. If real coupled modes are preferred, the use of tangential attachments modes is also an option but another technique is proposed. It consists in including in the basis the static response to a first-order approximation of the friction forces around the real coupled modes. Applied to the present disc/pads model, the two enrichment options provide errors about 10 times smaller than the initial basis. This is an interesting gain since errors on divergence rates are reduced to about 0.01%. Finally, it must be noted that the computational time gains obtained with these enriched reduced bases remain very high compared with iterative solvers. The corresponding cost/precision ratio is thus optimized.

Figure 1 :

 1 Figure 1: Local cartesian coordinate system at one particular point of the interface

Figure 2 :

 2 Figure 2: The disc/pads finite element model

Figure 3 :

 3 Figure 3: Free disc modes : (a) A0-5; (b) A1-3; (c) R0-3; and, (d) C0-1

Figure 4 :

 4 Figure 4: 'Sliding' pin modes: (a) F1; and, (b) T1

Figure 5 :

 5 Figure 5: Real coupled modes: (a) A04-, (b) A04+ and (c) and (d) two F1-type

Figure 6 :

 6 Figure 6: Stability charts for µ = 0.4 and: (a) C b = 0; or (b) C b 0

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Coalescence of complex disc modes: (a)-(c) A04-/A04+; and, (b)-(d) A05-/A05+; black dots: C b = 0; gray dots: C b 0

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: Performance of the free CMS basis for the frictionless problem for: (a) disc modes; (b) F1 pin modes; and, (c) T1 pin modes

Figure 14 :Figure 15 :

 1415 Figure 14: Performance of the free CMS bases for the stability study for: (a) disc modes; (b) unstable pins modes; and, (c) other pins modes

Figure 16 :Figure 17 :

 1617 Figure 16: Performance of the mixed CMS bases for the stability study for: (a) disc modes; (b) unstable pins modes; and, (c) other pins modes

Figure 18 :

 18 Figure 18: Global comparison of the performances of the enriched bases for the stability study

Table 1 :

 1 Characteristics of the 66 first components modes, real coupled modes and complex modes for µ = 0.

	3521	-0.12
	3525 to 3528	-
	3535	-0.51
	3537	+0.88
	3552 to 3563	-1.76 to -1.07

Table 2 :

 2 CMS bases for the frictionless problem

	Max.	n D m	n P m
	Frequency all methods free method fixed and mixed method
	6 kHz	30	36	36
	7.5 kHz	60	72	36
	10 kHz	80	120	60

Table 3 :

 3 Number of components normal modes taken into account in CMS

Table 4 :

 4 Enriched CMS bases for the stability study (t.a.m.: tangential attachment modes, r.a.m.: rubbing attachment modes)

Table 6 :

 6 Size and computational time gains of the four retained bases

		n r	G n t r [s] G t
	free + t.a.m.	428 72	142 23
	free with r.a.m. 200 155 108 30
	real + t.a.m.	294 105 124 26
	real + 1 st -ord. 132 234 113 29


S Besset

L Jézéquel Ltds

Équipe D2s

Optimization of Structural Dynamic Behaviour Based on Effective Modal Parameters

Keywords: Modal analysis, Optimization, Modal criteria

Optimization of complex structures often leads to high calculation costs. Indeed, the structure has to be frequently reanalysed in order to update the optimization criteriums. We propose an optimization method based on effective modal parameters. These parameters are close to the modal matrices used for the modal analysis of a structure. Thus, once the structure has been analysed, it becomes very easy to calculate optimization criteria. First, we will explain the modal analysis that we will use in this paper. A modal model will be used to analyse the hollow parts of the structure. The modal analysis of the whole structure will be performed using substructuring and "double modal synthesis" proposed by Jezequel. Secondly, we will explainthe how to obtain effective modal parameters and their use for optimization. Finally, we will show the efficiency of these parameters through the optimization of a complex structure, using two types of optimization methods.

Introduction

Using matrices resulting from the modal analysis of a structure in order to optimize the dynamic behaviour of this structure has multiple advantages. Once the structure has been analyzed, optimization criteria become very easy to compute. Secondly, it is possible to link the optimization of the structure and the modal matrices resulting from the modal analysis associated with the resonance phenomena. Thus, it is possible to identify the causes of the problems to solve. The criteria we will use in this paper have been developed by P. Lemerle [START_REF] Lemerle | Optimisation des structures selon des critères imposés par la discrétion acoustique des navires[END_REF], using the Craig & Bampton method for structure analysis. The doctoral dissertation of P. Lemerle is derivated in the present paper. Similar optimization methods have already been used by Suweca [START_REF] Suweca | Optimal structural design with damping constraint limitations[END_REF] in the case of structural designs with damping constraint limitations.

To optimize the structure, we will modify properties of the hollow parts of the structure, but other properties could be changed. The example we propose is a complex structure, including plates and stiffeners, that we will optimize through two types of methods. First, a global method, using a genetic algorithm. Secondly, a local optimization will be processed, using a method based on the gradient method and using Kühn & Tucker conditions.

Effective modal parameters

In order to compute the effective modal parameters, we will first explain the modal analysis of the structure used in this study.

Modal analysis of a structure

The structure used in this study is a complex structure including hollow parts and plates. It is made of formed steel constituting its skeleton, as shown in figure 1. Plates are attached to this skeleton. The geometry of the structure is similar to the geometry of a car, in order to demonstrate that the methods we propose can be used in an industrial context.

[Figure 1 about here.] Hollow parts of the structure -which could be stiffeners for example -are analyzed using a model we proposed in [START_REF] Besset | Modal analysis of hollow parts of a structure[END_REF]. These elements constitute the skeleton of the structure. This modelization leads to modal matrices that can be assembled like finite elements matrices. The main characteristic of this method is that it produces matrices including only generalized degrees of freedom. There may remain nodal degrees of freedom in order to assemble the hollow parts with other structures, but boundaries between the elements constituting the hollow part only comprise generalized degrees of freedom. It is possible to assemble these elements because of the choice of modes used for the modal analysis of the elements. We now summarize the method used in [START_REF] Besset | Modal analysis of hollow parts of a structure[END_REF]. Figures 2 and3 show the geometry of the hollow parts we consider. The hollow part shown in figure 2 is split into several elements shown in figure 3.

[Figure 2 When the modal analysis of one element is complete, two nodal degrees of freedom will remain, in order to couple the hollow part with another structure. The behavior of an element can be described with the equation of motion:

-ω 2 M LL M LR M RL M RR + K LL K LR K RL K RR u L u R = f L f R (1)
where mass and stiffness matrices M and K are split into left and right degrees of freedom L and R. Vectors u and f are split the same way. u L and u R are then split into degrees of freedom that will be expressed as a function of generalized degrees of freedom u g L and u g R , and the others u k L and u k R -marked "retained" on figure 3:

u L = u k L u g L , u R = u k R u g R (2)
Let φ R be the matrix of the modes of the element when one node of the right side is fixed, as shown in figure 4.

[Figure 4 about here.] Φ R is the modal matrix corresponding to right nodes of the element -it is a part of matrix φ R . According to the Craig & Bampton theory [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF], displacements of these right nodes can be expressed as follows:

u g R = Φ R q R + Ψ R u k R ( 3
)
where Ψ R is the matrix of the static modes, corresponding to the rigid body modes of the right side of the element.

Analogous matrices are defined for the left side of the element. Let φ L be the matrix of the modes of the element when one node of the left side is fixed. Φ L is the modal matrix corresponding to left nodes of the element. Displacements of these left nodes can be expressed as follows:

u g L = Φ L q L + Ψ L u k L (4)
Thus, displacements u can be expressed as a function of generalized degrees of freedom q L and q R :

       u k L u g L u k R u g R        =     I 0 0 0 Ψ L Φ L 0 0 0 0 I 0 0 0 Ψ R Φ R            u k L q L u k R q R        (5)
Notice that as left degrees of freedom do not depend on right modes, and right degrees of freedom do not depend on left modes, assembling modal elements will be possible.

In order to assemble two elements, let write the relations between the nodal degrees of freedom of the left and right side of two elements. Degrees of freedom of the right side of the first element, called u k R1 and u g R1 , must match with degrees of freedom of the left side of the second element, called u k L2 and u g L2 :

u k R1 = u k L2 u g R1 = u g L2 (6)
Equation 6 can be written for elements n and n + 1:

u k Rn = u k Ln+1 u g Rn = u g Ln+1 (7)
Equations 6 leads to the expression of degrees of freedom of the left part of the second element as a function of degrees of freedom of the first element:

u k L2 = u k R1 q L2 = ΦL Φ R q R1 + (Ψ R -Ψ L ) u k R1 (8)
which can be written for elements n and n + 1:

u k Ln+1 = u k Rn q Ln+1 = ΦL Φ R q Rn + (Ψ R -Ψ L ) u k Rn ( 9
)
where ΦL is a pseudo-inverse matrix of Φ L :

ΦL = Φ T L Φ L -1 Φ T L ( 10
)
Equations 8 allow the assembling of N elements, using transfer matrices T n :

K tot = N n=1 T T n KT n (11)
M tot = N n=1 T T n MT n (12)
Plates are assembled with the skeleton of the structure through the nodal degrees of freedom remaining from the modal analysis of the hollow parts of the structure. In this paper, only one plate will be used, located on the top of the structure. The entire structure is then analyzed using the "double modal synthesis" method proposed by Jezequel [START_REF] Jezequel | Component modal synthesis methods based on hybrid models, part i : Theory of hybrid models and modal truncation methods[END_REF][START_REF] Jezequel | Component modal synthesis methods based on hybrid models, part ii : Numerical tests and experimental identification of hybrid model[END_REF]. This method uses "branch modes" to describe the behaviour of the boundaries between substructures. In this paper, will used these "branch modes" to discribe the behaviour of the skeleton of the structure.

Thus, mass and stiffness matrices of the structure can be split into degrees of freedom concerning plates and degrees of freedom concerning hollow parts of the structure. Degrees of freedom concerning hollow parts include generalized degrees of freedom resulting from the modal analysis of the elements constituting the skeleton, that will be denoted q Hc , and generalized degrees of freedom resulting from the double modal synthesis, which will be denoted q Hb . We denote q H = q Hc q Hb . Vector q P corresponds to the generalized degrees of freedom of the plates. Hence the motion equation:

-

ω 2 M HH M HP M P H M P P + K HH K HP K P H K P P q H q P = f H f P (13)
Generalized degrees of freedom used in the motion equation 13 are linked to nodal degrees of freedom through the following equations:

   u P = Φ P q P + Ψ P u Hb u Hc = Φ Hc q Hc + Ψ Hc u Hb u Hb = Φ Hb q Hb (14)
where u P are the nodal displacements of the plates, u Hc are the nodal displacements of the hollow parts corresponding to the generalized degrees of freedom q Hc , and u Hb are the nodal displacements on the boundaries, associated to the branch modes.

Matrix Φ P is the modal matrix of the fixed modes of the plates. Ψ P is the matrix of the static modes of the plates, as in the Craig & Bampton method [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF].

Matrices Φ Hc and Ψ Hc are modal matrices resulting from the analysis of the hollow part described in [START_REF] Besset | Modal analysis of hollow parts of a structure[END_REF].

Matrix Φ Hb is the matrix of the "branch" modes of the structure [START_REF] Jezequel | Component modal synthesis methods based on hybrid models, part i : Theory of hybrid models and modal truncation methods[END_REF][START_REF] Jezequel | Component modal synthesis methods based on hybrid models, part ii : Numerical tests and experimental identification of hybrid model[END_REF].

Obtaining effective modal parameters

Effective modal parameters used by P. Lemerle [START_REF] Lemerle | Optimisation des structures selon des critères imposés par la discrétion acoustique des navires[END_REF] are parameters that link degrees of freedom that are submitted to a displacement -excited degrees of freedom -and degrees of freedom whose displacements are to be minimized.

In order to obtain the effective modal parameters, we will use the modal analysis method proposed in section 2.1, but we will separate the degrees of freedom that will be excited -these degrees of freedom will remain nodal degrees of freedom.

In this paper, the excitations will concern some of the degrees of freedom we denoted u Hb . These degrees of freedom will be called u He . Hence equation 15, coming from equation 14:

   u P = Φ P q P + Ψ P u Hb + Ψ P e u He u Hc = Φ Hc q Hc + Ψ Hc u Hb + Ψ He u He u Hb = Φ Hb q Hb (15) Thanks to the orthogonal properties of the modes used in the modal analysis, equation 13 can be written as follows, considering the damping matrix which we assume to be diagonal. Indeed, we use the modal damping assumption (also called Basile's hypothesis in French terminology), hence the following equation:

    -ω 2     M EE M EHc M EHb M EP M HcE M HcHc M HcHb M HcP M HbE M HbHc m Hbk M HbP M P E M P Hc M P Hb m P k     + iω     0 0 0 0 0 c Hck 0 0 0 0 c Hbk 0 0 0 0 c P k     +     K EE K EHc K EHb 0 K HcE K HcHc K HcHb 0 K HbE K HbHc k Hbk 0 0 0 0 k P k                u He q Hc q Hb q P        =        f E f Hc f Hb f P        (16)
where matrices

[m Hbk ], [m P k ], [k Hbk ], [k P k ], [c Hbk ], [c P k ] and [c Hck ] are diago- nal matrices.
The analysis of the role played by the structural damping is not the purpose of the paper, that is why we chose a very simple way to take it into account.

To obtain the effective modal parameters, u P must be expressed as a function of f P and u He . Let us express one of the last lines of equation 16:

-ω 2 M k P E u He + M k P Hc q Hc + M k P Hb q Hb + m P k q k P + iωc P k q k P + k P k q k P = f k P (17)
M k P E , M k P Hc , M k P Hb are the k th lines of matrices M P E , M P Hc , M P Hb . f P can be expressed in function of f P as follow:

f P = Φ T P f P (18)
Equation 17 becomes:

q k P = Φ kT P f P + ω 2 M k P E u He + M k P Hc q Hc + M k P Hb q Hb -ω 2 m P k + iωc P k + k P k (19)
where Φ k P is the k th column of Φ P . Equation 19 becomes:

u P = k Φ k P q P k + Ψ P u Hb + Ψ P e u He = k Φ k P Φ kT P -ω 2 m P k + iωc P k + k P k f P +   k   ω 2 Φ k P M k P E -Mk P Hc Ψ He -ω 2 m P k + iωc P k + k P k   + Ψ P e   u He +   k   ω 2 Φ k P Mk P Hb -Mk P Hc Ψ Hc -ω 2 m P k + iωc P k + k P k   + Ψ P   u Hb + k -ω 2 Φ k P Mk P Hc -ω 2 m P k + iωc P k + k P k u Hc (20)
where M k P Hb = Mk P Hb Φ Hb and M k P Hc = Mk P Hc Φ Hc . Matrix ΦHc is a pseudoinverse matrix of Φ Hc .

Two effective modal parameters can be deduced from equation 20. First, the dynamic flexibility matrix G is given by:

G(ω) = k Φ k P Φ kT P -ω 2 m P k + iωc P k + k P k (21)
It corresponds to the relation between a force applied on a plate and the displacements it causes.

Secondly, ignoring the static terms corresponding to the boundaries, the transmissibility matrix T is given by:

T(ω) = k   ω 2 Φ k P M k P E -Mk P Hc Ψ He -ω 2 m P k + iωc P k + k P k   (22)
It is possible to take the static terms into account, but we are more interested in the modal sum k (), wich allows to analyze the contribution of each mode of the structure, hence the choice of the transmissibility matrix of equation 22.

This equation is linked to the relation between the excitation and the displacements it causes on the plate. Note that the excitation is a displacement excitation. They can be rewritten using effective parameters G and T as follows:

G(ω) = k 1 1 -ω ω k 2 + 2iξ k ω ω k

Gk

(23)

T(ω) = k ω ω k 2 1 -ω ω k 2 + 2iξ k ω ω k Tk ( 24
)
where

Gk = Φ k P Φ kT P ω 2 k m P k (25)

Optimization based on effective modal parameters

Criteria used for optimization

In this section, we will deduce criteria from the flexibility and transmissibility matrices proposed in section 2.2. The sums k () that appear in these matrices correspond to a superposition of modes. Thus, optimization criteria can be written as follow:

C G = max k Φ k P Φ kT P ω 2 k m P k (27)
C T = max k Φ k P M k P E -Mk P Hc Ψ He m P k (28)
where the norm |x| is the maximal component of matrix x. Considering these criteria, it is possible to optimize the structure. Moreover, obtaining the value k max allows us to understand which mode is responsible for the value of the criteria.

Methods used for optimization

We will show examples using the criteria we propose through a two level strategy for optimization. First, a genetic algorithm will be used. Secondly, a local optimization will be processed in order to obtain more precise results.

Genetic optimization

Genetic optimization is often used in the case of multi-objective optimization problems. The aim of this section is to generate a Pareto front in order to obtain various optimized solutions for the optimization problem. Overviews of multi-objective optimization genetic algorithms can be found in [START_REF] Stadler | A survey of multicriteria optimization, or the vector maximum problem[END_REF][START_REF] Fonseca | An overview of evolutionary algorithms in multiobjective optimization[END_REF][START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learneing[END_REF]. The genetic algorithm that will be used in this section is quite simple. Let α i be the p parameters to optimize. A 1 is a first matrix of N sets of parameters:

A 1 =        α 1 1 α 1 2 • • • α 1 p-1 α 1 p α 2 1 α 2 2 • • • α 2 p-1 α 2 p . . . . . . . . . . . . . . . α N -1 1 α N -1 2 • • • α N -1 p-1 α N -1 p α N 1 α N 2 • • • α N p-1 α N p        (29)
The lines of matrix A 1 are the N sets of parameters chosen for the first iteration of the algorithm. The values of the criteria we propose are computed using these sets of parameters. The mass of the structure can also be a third criterion, for the structure mustn't become too heavy. A set of α i is then selected if no other set gives better results on both criteria -and the mass. The other sets of α i are combined in order to create another matrix A 2 of N sets of parameters that will include the selected sets. Figure 5 sums up the algorithm used for this optimization.

[Figure 5

about here.]

Once the optimal sets of α i have been found, the results can be plotted on a pareto diagram, which will be done bellow. Various "optimal solutions" are found, and the local optimization can be made from one of these "optimal solutions".

Local optimization

The genetic algorithm proposed in section 3.2.1 gives various results, but the best value for each α i is not precisely determined. In order to obtain more accurate results, we propose a local method.

Three criteria are to be considered for the optimization. C G and C T have been defined in section 2.2. The third criterion will be the mass of the structure, denoted C M .

The optimization problem can then be written:

Minimize C M (α) assuming CG (α) ≤ 0 CT (α) ≤ 0 (30)
where α = [α 1 , α 2 , . . . , α p ]. CG and CT are directly linked with C G and C T . The method used in this paper requires that the criteria be derivated, which is not possible with C G and C T . Such criteria are defined as follow:

CG = 1 4 log k Φ k P Φ kT P -ω 2 k m P k 4 -c g (31) CT = 1 4 log k   Φ k P M k P E -Mk P Hc Ψ He m P k   4 -c t ( 32
)
where c g and c t are the objectives of the criteria. CG and CT have almost the same minima and extrema as C G and C T . To prove this, let us consider the function f p as follows:

f p : x -→ p a 1 (x) p + a 2 (x) p + • • • + a n (x) p (33)
where a q (x) > 0 ∀q ≤ n ∀x. It is possible to write:

lim p→+∞ f p (x) = max [a 1 (x), a 2 (x), . . . , a n (x)] (34)
Indeed:

p max [a 1 (x), a 2 (x), . . . , a n (x)] p < f p (x) < p n max [a 1 (x), a 2 (x), . . . , a n (x)] p ∀p ∈ N. ( 35
)
Equation 35 leads to:

max [a 1 (x), a 2 (x), . . . , a n (x)] < f p (x) < p √ n max [a 1 (x), a 2 (x), . . . , a n (x)] ∀p ∈ N. ( 36
)
As lim

p→+∞ p √ n = 1, we can deduce lim p→+∞ f p (x) = max [a 1 (x), a 2 (x), . . . , a n (x)].
Thus, if p is high enough, CG and CT have almost the same minima and extrema as C G and C T . In this section, we will use CG and CT instead of C G and C T , for these criteria may be derivated, whereas C G and C T cannot.

The optimization method we use in this paper is based on the Kühn & Tucker conditions, which can be written as follow:

∃ (λ 1 (α), λ 2 (α)) , λ 1 (α)∇ CG (α) + λ 2 (α)∇ CT (α) + ∇C M (α) = 0 and λ 1 (α) ≥ 0 λ 2 (α) ≥ 0 (37)
The problem we propose requires equation 37 to be written as follows:

∀i = 1, 2, 3, . . . , p λ 1 (α) ∂ CG ∂αi (α) ∂CM ∂αi (α) + λ 2 (α) ∂ CT ∂αi (α) ∂CM ∂αi (α) = -1 (38)
In this equation, we assume ∂C M ∂α i (α) = 0. A relaxation parameter γ is then introduced. Equation 38 becomes:

α i = γα i -(1 -γ) λ 1 (α) ∂ CG ∂αi (α) ∂CM ∂αi (α) + λ 2 (α) ∂ CT ∂αi (α) ∂CM ∂αi (α) α i ( 39
)
This equation leads to the following recurrence equation:

α k+1 i = γα k i -   (1 -γ)   λ 1 (α) ∂ CG ∂αi (α) ∂CM ∂αi (α) αi=α k i + λ 2 (α) ∂ CT ∂αi (α) ∂CM ∂αi (α) αi=α k i     α k i ( 40
) This method is similar to the gradient method, but the descent direction here is:

d k = -(1 -γ)R i (α)α k i (41)
where:

R i (α) = 1 + λ 1 (α) ∂ CG ∂αi (α) ∂CM ∂αi (α) αi=α k i + λ 2 (α) ∂ CT ∂αi (α) ∂CM ∂αi (α) αi=α k i (42)
The convergence condition can then be written as for the gradient method. R : α -→ {R i (α)α i } must satisfy the following equations:

∃κ > 0 ∀u, v R(v) -R(u), v -u ≥ κ ||v -u|| 2 (43) ∃ξ > 0 ∀u, v R(v) -R(u) ≤ ξ ||v -u|| (44)
It is then possible to find a value for γ. The algorithm converges if the following relation is satisfied:

1 - 2κ ξ 2 < γ < 1 (45)
γ must be chosen close to 1 if we want to be sure that the algorithm converges. The smaller γ is, the faster the algorithm will converge, but it may not converge at all if γ is too small. . . For each step of the algorithm, λ i (α) is being computed. The method used to obtain these λ i (α) is now explained. Let us take ∆α i , such as:

CG (α 1 + ∆α 1 , α 2 + ∆α 2 , . . . , α p + ∆α p ) = 0 CT (α 1 + ∆α 1 , α 2 + ∆α 2 , . . . , α p + ∆α p ) = 0 (46)
This equation allows us to express ∆ CG (α 1 , α 2 , . . . , α p ) as follows:

∆ CG (α 1 , α 2 , . . . , α p ) = CG (α 1 + ∆α 1 , α 2 + ∆α 2 , . . . , α p + ∆α p ) -CG (α 1 , α 2 , . . . , α p ) = -CG (α 1 , α 2 , . . . , α p ) = p i=1 ∂ CG ∂α i ∆α i (47)
Replacing ∆α i in equation 47, CG can be written as follows:

CG (α 1 , α 2 , . . . , α p ) = (1 -γ) p i=1 ∂ CG ∂α i 1 + λ 1 ∂ CG ∂αi ∂CM ∂αi + λ 2 ∂ CT ∂αi ∂CM ∂αi α k i (48)
In order to simplify the notations, we denote

λ 1 = λ 1 (α), λ 1 = λ 1 (α) and ∂ Cj ∂α i = ∂ Cj ∂α i (α) αi=α k i .
The same equation can be written for CT , which allows us to write the following equation:

        p i=1 ∂ CG ∂αi 2 ∂CM ∂αi α k i p i=1 ∂ CG ∂α i ∂ CT ∂αi ∂CM ∂αi α k i p i=1 ∂ CT ∂α i ∂ CG ∂αi ∂CM ∂αi α k i p i=1 ∂ CT ∂αi 2 ∂CM ∂αi α k i         λ 1 λ 2 =            CG 1 -γ - p i=1 ∂ CG ∂α i α i CT 1 -γ - p i=1 ∂ CT ∂α i α i           
(49) This equation allows us to obtain (λ 1 , λ 2 ) for each step of the optimization.

Results

The method proposed in this paper has been tested on a complex structure including hollow parts and plates. The parameters we choose to optimize are linked to the geometry of the hollow parts. Indeed, we will optimize D and λ, as shown in figure 6.

[Figure 6 about here.]

The optimization methods we will use in the next sections will use D and λ as parameters. The hollow parts of the structure shown in figure 1 are split into 8 parts, and each part is optimized with optimal values of D and λ. Thus, 16 parameters are to be optimized.

During the optimization, matrices M and K are updated in order to take into account the evolution of the geometry of the hollow parts of the structure. Actually, modes concerning these hollow parts need to be re-computed for each iteration. However, only one element for each set of parameters has to be re-computed, thanks to the use of our substructuring method to describe the hollow parts of the structure. Parametric methods can also be used to solve numerical difficulties.

Structure used for optimization

The structure used for the optimization is given in figure 1. It is made of hollow parts and plates, that have been analyzed using the methods proposed in section 2.1.

Because of the calculation costs, we chose a structure which is not too large, i.e. 2.5 meters long, as given in figure 1.

Analysis of the criteria

The analysis of criteria C G and C T shows which modes are responsible for the displacements needing to be reduced. Figure 8 shows the values of C k G , which are part of criterion C G :

C k G = Φ k P Φ kT P ω 2 k m P k (50)
The same analysis can be made for C k T :

C k T = Φ k P M k P E -Mk P Hc Ψ He m P k (51)
Figure 7 shows the values of C k T , which are part of criterion C T .

[Figure 7 about here.]

[Figure 8 about here.]

Figures 7 and8 show that criteria C G and C T do not necessarily depend on the same modes. In this example, we only analyze the criteria for the first 50 modes of the structure. Figure 7 shows that the 45 th mode is mostly responsible for the value of criterion C G , whereas figure 8 shows that the 6 th , 7 th , 22 th and 23 th modes are responsible for the value of criterion C T . Thus, it is necessary to take these two criteria into account to optimize the structure.

Results of genetic optimization

In this section, we present the results obtained using the genetic algorithm given in figure 5. Figures 9,10,11 and 12 show the pareto diagrams of the sets A i . Units on x and y axis are not important because it depends on the values of c g and c t . The excitation points are located on the hollow parts, whereas the displacements to be reduced are located on the plate on the top of the structure.

[Figure 9 These results show that the algorithm is able to produce a set of optimal parameters. These optimal parameters form a curve including points that are the closest to the origin of the graphs. The pareto diagram is shown in 3D, but because it is not easily understood, we have also shown three 2D diagrams.

Figure 13 shows the evolution of the number of pareto points as a function of the progress of the optimization. At the end of the algorithm, all points of the set A N are "pareto points" -i.e. optimal points.

[Figure 13 about here.]

Results of local optimization

A local optimization is now processed, using the same structure and the same excitation points. Figure 14 shows the square speed of a point of the plate, whereas figure 15 shows the evolution of the criteria during the optimization.

[Figure 14 In figure 15, the ratio between the criteria and their initial values is plotted. That is why the first values are 1 in the three cases. Criteria do not change much because the values chosen for the parameters α i are close to the optimal values. This is a condition for the algorithm to work well. Notice that the algorithm is not always able to converge because some conditions on the parameters have been set. The parameters to be optimized must remain in an interval of acceptable values.

In figure 14, the optimized value for the square speed is not much better than the initial value because the parameters have not changed much. Nevertheless, this value is not bad and the mass of the structure has decreased, which is why the optimal solution is better than the initial one.

Conclusion

The optimization method we propose in this paper is based on criteria that allow us to analyze the origin of the problems we want to solve -reducing noise in complex structures. These criteria have been used with two kinds of optimization algorithms and the results prove that they are quite efficient and able to represent th vibration level of the structure.

Examples given in this paper are based on an excitation located on a point of a hollow part. It is possible to do the same analysis based on another excitation. Criteria will then be different, but the results should be quite good. The algorithms proposed here are quite easy to use with the structures we chose.

  about here.] [Figure 3 about here.]

  c P k = 2ξ k √ k P k m P k and ω k = k P k m P k . Matrices Gk and Tk are called effective modal parameters.

  about here.] [Figure 10 about here.] [Figure 11 about here.] [Figure 12 about here.]

  about here.] [Figure 15 about here.]

Figure 1 :Figure 3 :Figure 4 :

 134 Figure 1: Structure to be optimize

Figure 5 :Figure 6 :Figure 14 :Figure 15 :

 561415 Figure 5: Algorithm used for genetic optimization

List of Figures