Thomas Philippe
email: thomas.philippe@etu.univ-rouen.fr

Didier Blavette

Didier Gpm Blavette

; Gpm

Nucleation pathway in coherent precipitation

Keywords: solid-state transformation, precipitation, nucleation minimum energy pathway solid state phase transformations, precipitation, non-classical nucleation, minimum energy pathway

published or not. The documents may come

Introduction

Numerous solid state transformations among them precipitation proceed via first-order phase transitions. According to the supersaturation, the early stages of precipitation in alloys can be classically accounted for either by nucleation and growth of embryo or spinodal decomposition for which no barrier is involved. The early developments of the thermodynamics of first-order phase transitions were performed by Gibbs [1] and van der Waals [START_REF] Van Der Waals | Lehrbuch der Thermodynamik[END_REF][START_REF] Rowlinson | Translation of J.D. van der Waals[END_REF]. The first formulation of the nucleation-growth model was done by Gibbs [START_REF] Volmer | Kinetik der Phasenbildung[END_REF][START_REF] Ya | Kinetic Theory of Liquids[END_REF][START_REF]Nucleation[END_REF][START_REF] Zettlemoyer | [END_REF][START_REF] Gunton | Phase Transitions and Critical Phenomena[END_REF][START_REF] Skripov | Metastable Liquids[END_REF]. Cahn and Hilliard re-expressed the van der Waals method to develop a continuum theory of nucleation [START_REF] Cahn | [END_REF]11]. Many developments and refinements were brought about to the early theories of pioneers like Becker [12]. Let us mention the models of Zeldovich and Turnbull which provide us with a latent nucleation time [13,14]. When the equilibrium phase has a different structure as the parent solid solution, the system often prefers to form a transient coherent phase that has the same structure as the initial solid solution (coherent nucleation). This is well known in Al base alloys, notably in AlCu where coherent GP zones nucleate prior the precipitation of the equilibrium phase Al 2 Cu [15]. It is well accepted that these GP zones are almost pure in copper. However, in other systems, and outside the spinodal region, dilute and diffuse clusters are observed to form before the equilibrium phase precipitates. Such nonclassical nucleation processes were observed both in metals (Fe-Cu, Fe-Cr) and in semiconductors (Si-B, Si-As) using atom probe tomography (APT) [16,17]. Because of its ultimate spatial resolution, APT has played a major role in the investigation of the early stages of decomposition [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. Although the system is in a metastable state (e.g. outside the spinodal regime of decomposition where the system is in an unstable state), nuclei were observed with a lower solute content than the equilibrium phase and found to enrich with time before reaching the expected composition. Interfaces sometimes appear diffuse. Diffusion as well as thermodynamics control kinetics pathway during precipitation [32][33][34]. In this paper, diffusion effects have been ignored. Only arguments based on thermodynamics are taken into account.

The nucleation of embryos that have a low solute content may appear surprising when considering the only contribution of the driving force for nucleation ∆f. This force is clearly shown in figure 1 to decrease when the solute content in the nucleus is smaller. However, this can be counterbalanced by the lower interfacial energy (σ) that occur between diffuse nuclei and the parent phase. σ is indeed proportional to the concentration gradient (∇c) 2 is consequently much lower for dilute and diffuse nuclei. Dilute nuclei also leads to a smaller mismatch between phases that reduces the misfit elastic energy that also acts against nucleation. In other terms, in spite of the lower driving force, the barrier of nucleation of such dilute and diffuse solute-enriched zones can be lower than that of the expected coherent phase. This article aims at predicting the optimum composition of nuclei that minimises the nucleation barrier. Kinetics arguments related to diffusion were not considered in this approach. Another question to which this article wishes to answer is also how the composition of nuclei evolves towards equilibrium with increasing size.

The diffuse interface theory of non-classical nucleation developed by Cahn and

Hilliard has been the first to deal with this issue [START_REF] Cahn | [END_REF]11]. The theory is based on the determination of the radial concentration profile of the critical nucleus minimizing the nucleation barrier. They found that the classical nucleation theory was only operative for small supersaturation. Composition profile of a critical nucleus as a function of the nominal composition (c 0 ) from the equilibrium composition (c α ) to the spinodal one (c s ) was investigated in regular solutions (T/Tc=0.8656). They found that the nucleus composition (at the core of the nucleus) approaches the matrix composition near the spinodal line. The critical radius of nucleation was found to diverge to infinity not only near c α but also close to the spinodal line c s . The strong point of this diffuse interface theory is that no hypothesis about the chemical homogeneity of nuclei is required. More recently, Schmelzer et al. developed a simpler model (the generalized Gibbs theory) based on the same ideas as Cahn and Hilliard, i.e the minimization of the nucleation barrier [35]. In contrast to Cahn, the authors considered nuclei with an homogeneous In this work, we have focused on the nucleation of precipitates in solids. We have therefore extended the generalized Gibbs model developed by Schmelzer et al. [35,[START_REF] Schmelzer | On the Dynamics of First-Order Phase Transitions in Multi-Component Systems[END_REF] so that to include the contribution of the misfit elastic energy [37]. The latters have studied the influence of elastic stress on crystallisation in glasses [38]. Here, we have focused on the influence of misfit on solid state precipitation. We have confronted predictions to experiment results provided by atom-probe tomography. The latter is one of the very rare instruments able to measure the composition of small embryo (1 nm) in the early stages of nucleation. The model is only based on thermodynamics arguments.

The influence of diffusion and coupling effects of diffusion fluxes on the kinetics pathway was not considered. Nucleation was modelled as a function of the composition and radius of nucleus. The composition of embryo was assumed to be uniform and was not constrained to be that of the equilibrium phase contrary to the classical nucleation theory. The chemical composition (c*) of the critical nucleus as a function of the nominal composition was computed through the determination of the saddle point of the nucleation barrier W(R,c) in the R,c space. A regular solution was considered in order to compute the driving force for nucleation. Isotropic elastic energy was considered (Eshelby theory of inclusions [37]). The most probable trajectory, the minimum energy pathway, along the free energy surface W(R,c) was computed numerically. This made it possible to predict how the composition evolves as a function of the embryo size in the early stages of nucleation. Confrontation of predictions to experiments dealing with the early stages of precipitation in Si-B, Fe-Cr and Fe-Cu systems has been carried out. Interfacial and elastic energies have been assumed to be isotropic. The nucleus composition has been assumed to be uniform. The nucleation of the new phase called hereafter β has been investigated as a function of the size (spherical nucleus of radius R) and composition (nucleus composition c). As for classical nucleation theory, the composition of the parent solid solution was assumed to be nearly constant during nucleation (negligible volume fraction of β phase). The free energy change due to the formation of a nucleus can be classically expressed as:

( ) F f V S γ σ ∆ = -∆ - + (1)
f ∆ (in J/m 3 ) is the driving force for nucleation (Fig. 1) is defined as following:

0 0 0 ( ) ( ) ( ) c c f f f c c c f c c = ∂ ∆ = + - - ∂ ( 2
)
where c is the nucleus composition (i.e. atomic fraction of B), c 0 the nominal one and the function f is the free energy (fig. 1). For the sake of simplicity, an isotropic inclusion has been considered. The elastic energy γ, as given by Eshelby [37] can be written as proportional to the square function of the lattice misfit between the nucleus β and the parent phase α. The misfit being proportional to the composition difference between both phases (Vegard law), one can simply express γ as: The interfacial energy was written as proportional to the concentration gradient squared:

2 0 ( ) c c γ γ = - % (3)
2 0 ( ) c c σ σ = - % (4)
with σ% =ε/d 2 , ε the pair interaction energy and d the interface thickness. This latter parameter (d) may account for the diffuseness of interfaces.

Nuclei were assumed to be spherical (radius R) so that the change of free energy reads:

( )

3 2 2 2 0 0 4 ( ) 4 ( ) 3 F R f c c R c c π γ π σ ∆ = - ∆ - - + - % % (5)
The free energy surface, ∆F(R,c), for a given nominal composition (c 0 ) acts as a thermodynamic barrier over which the growing nucleus must pass through before growth. The major nucleation flux generally passes through the saddle point [12]. The location of this saddle point is defined by the two following equations:

2 0 2 0 2 ( ) 0 ( ) c c F R R f c c σ γ - ∂∆ = ⇔ = ∂ ∆ - - % % ( 6
) 0 0 0 2 ( ) 2 ( ) 0 3 F R f c c c c c c γ σ ∂∆ ∂∆   = ⇔ - - - + - =   ∂ ∂   % % (7)
The composition of the critical nucleus (c*) can be derived by combining equations ( 6) and ( 7):

8 0 * 0 3 ( * ) * c c f f c c c c c γ = ∂∆ ∆ + -= ∂ - % (8)
The expression of ∆f (section b) makes it possible to find c*. Once the critical composition (c*) is determined, the critical radius (R*) can be derived from equation [START_REF]Nucleation[END_REF]. It is worth mentioning that equation ( 8) does not depend on the interfacial energy factorσ% . As a result, surprisingly, c* is independent of σ% (eq.8). The work required to form a nucleus (nucleation barrier W) is defined as the free energy change associated with the formation of a critical nucleus (R*,c*):

( )

3 6 0 0 2 2 0 ( * ) 16 ( *, ) 3 ( * ) c c W F c c f c c σ π γ - = ∆ = ∆ - - % % (9)
As the interfacial term of W grows as the 6 th power of the nucleus composition, it may be anticipated that the nucleation of nuclei of low solute content (smaller W) may be favoured although the driving force is smaller for such dilute nuclei. The model will be applied in the next section to a regular solution for the sake of simplicity. More realistic functional can however be used.

b. Non-classical nucleation in a regular solution

Let us consider a regular solution subject to unmixing. The free energy (f) can be expressed as:

2 (1 ) ln (1 ) ln(1 ) B f c c c c c c k T η = -+ + - - ( 10
)
where c is the atomic fraction of solute atoms and η the reduced temperature defined as T/T c with T the temperature and T c the critical temperature. The critical temperature T c is related to the interaction parameter ε by The binodal line (df/dc=0 -solubility limits c α ) and the spinodal line (c s ,

d 2 f/dc 2 =0
) in the (η,c) phase diagram are given respectively by:

1 2 ln (1 2 ) c c c α α α η   - = -     (12) (1 ) 4 s s c c η
- = (13)
Both curves classically coincide at the critical point (η=1, c=0.5). As an illustration, let us choose a system for a reduced temperature η=0.3 (solubility limit c α ~0.0013). The spinodal curve gives c s ~0.0816. We have focused on the range c α <c 0 <c s where precipitation proceeds via nucleation (spinodal decomposition has not been considered). The driving force for the formation of a nucleus with a given composition c (∆f(c), Fig. 1) can be computed as a function of the nominal composition (c 0 ).

Combining equations ( 2) and ( 10) leads to:

0 . ( , ) B f k T c c ρ ∆ = Γ (14) with 2 0 0 0 0 2 1 ( , ) ( ) ln (1 ) ln 1 c c c c c c c c c c η   - Γ = - - + -   -   ( 15
)
where ρ is the atomic density (at/m 3 ). Note that the driving force does not depend on the free energy of pure A, B components (F A and F B ). In the classical nucleation theory, Here, the composition of the nucleus is not constrained to be that of the equilibrium coherent phase (c β ). The consequence is that we do not consider any local equilibrium at the α/β interface at stated in classical growth theories. We consider that the nucleus composition is the one that minimises the nucleation barrier. However, all values of c* are not valid. Indeed, a phase separation may only occur if the thermodynamic driving force of the transformation ∆f is strickly positive (fig. 1). The minimum value (c m , fig. 1) of c* is given by ∆f=0, i.e. Γ(c,c 0 ) = 0 (eq.14 and 15). The latter equation has no analytical solution. For a given initial composition, c m increases when T/Tc increases. Once the composition of the critical nucleus c* is determined (equation ( 8)), the expressions of W and R * can be derived for c=c* from equation ( 6) and ( 9).

Properties of the critical clusters as predicted by this model are shown and compared to predictions of the classical nucleation theory in Fig. 2. We will first examine the situation where the elastic contribution is null (γ=0 -null misfit). The composition of the critical cluster is observed to drastically decrease with the initial solute concentration c 0 (Fig. 2a). For small supersaturations (c 0 close to c α ), c* is close to equilibrium (c β ~ 1) as stated in the classical theory. It decreases with increasing values of c 0 and tends towards c 0 when approaching the spinodal line (c 0 ~cs ). The driving force for nucleation (equation ( 14)) tends to zero when the spinodal curve is approached (fig. When the nominal composition approaches the spinodal line, c*→c s , Γ→0, R*→∞ and W→0. This shows that the transition between the nucleation-growth and the spinodal regimes is smooth in contrast to the predictions of classical nucleation theory.

It should be kept in mind that nucleation in solids may only occur if the driving force is larger than the elastic energy so that to ensure that the critical radius remains positive. In the classical case, this condition leads to an additional constraint on the nominal composition c 0 that should be larger than a minimum value given by:

2 0 0 ( , ) ( ) 0 B k T c c c c β β ρ γ Γ - - > % (16)
The properties of critical clusters are represented as a function of the initial solute concentration in Fig. 2 for different value of the elastic constant k (

/ B k k T γ ρ = % ).
k=1,2 have been used to study the influence of relatively high misfit on the properties of critical clusters. Typically, k=1 corresponds to a misfit of ~5.6% and k=2 to a misfit of ~7.9% (E=130 GPa, ρ=5*10 22 at/cm 3 and T=873 K). Again, for small supersaturation (c 0 small), c* tends to equilibrium (c β ). The minimum composition c 0 required for nucleation to start (positive effective driving force, equation ( 16)) is thus equivalent to the classical one. This composition corresponds to the solubility limit as given by the coherent bimodal composition (close to 0.01 for k=2, i.e., the coherent T/T c is equal to 0.43). As shown in Fig. 2a, the composition of the critical cluster increases drastically when the elastic energy parameter k increases. In other terms, the minimization of the 2b). The work required to form a critical nucleus (W) also increases with the elastic energy. It was also find that W vanishes and R* diverges when the initial solute concentration approaches the coherent spinodal line (Fig. 2c and2d). The coherent spinodal defined as d 2 f/dc 2 = -2γ% is located inside the incoherent spinodal line (d 2 f/dc 2 = 0), i.e. the spinodal limit is found for higher concentrations (minimum of W in Fig. 2c). As already demonstrated by Cahn and Hilliard, the classical nucleation theory tends to overestimate the work of formation of the critical nucleus. This difference could have a great influence on the calculation of the number density of nuclei (/m 3 ) as it depends on the exponential function of W.

Minimum energy pathway of nucleation

We propose here to determine the properties of the sub-and supercritical clusters. For that purpose, we consider that the evolution of the cluster properties (R,c) proceed along a valley of the thermodynamic potential (Fig. 3, i.e the 2D nucleation barrier). At the saddle point of the free energy surface W(R, c), the direction of negative curvature is related to the unstable variable (here R), while the direction of positive curvature is associated to the stable one (c). Whereas for many problems of physics the behaviour of stable variable is more interesting, it is not the case for nucleation. The minimum energy path is defined as the evolution of the cluster along the nucleation barrier following the path of the steepest descent from the saddle point. Numerical approaches to derive the minimum energy pathway of nucleation in multidimensional space were recently proposed (the nudged elastic band method [START_REF] Henkelman | Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points[END_REF][START_REF] Henkelman | A climbing image nudged elastic band method for finding saddle points and minimum energy paths[END_REF] and the string method [START_REF] Ren | Simplified and improved string method for computing the minimum energy paths in barrier-crossing events[END_REF][START_REF] Ren | String method for the study of rare events[END_REF]). Zhang et al. have investigated the computed the MEP of nucleation using the string method from the phase field equation [43]. They focused their work on the order parameter and the morphology of critical and equilibrium nuclei [START_REF] Zhang | Morphology of critical nuclei in solid state phase transformations[END_REF][START_REF] Zhang | Diffuse-interface description of strain-dominated morphology of critical nuclei in phase transformations[END_REF][START_REF] Zhang | Mathematical and Numerical Aspects of Phasefield Approach to Critical Morphology in Solids[END_REF].

Let us now consider the dimensionless reduced energy landscape Φ(r,c) defined as :

( , ) ( , ) F r c r c W ∆ Φ = ( 17
)
where W is the height of the nucleation barrier at the saddle point. Let us also introduce the reduced radius r = R/R * . For a given composition of the system (c 0 ), Φ defines the landscape of the energy barrier (fig. 3). Using equations (1) ( 5), ( 6), ( 9) and ( 17), the function Φ can finally be written as:

2 2 3 2 0 0 0 2 2 0 0 0 ( , ) ( ) ( ) ( , ) 2 3 ( *, ) ( * ) ( * ) f c c c c c c r c r r f c c c c c c γ γ ∆ - - - Φ = - + ∆ - - - % % (18)
Surprisingly, Φ(r,c) does not depend explicitly on the interfacial constant. There is however an implicit dependence of Φ as a function of R as r=R/R* and R* is controlled by σ. It is worth mentioning that there is no analytical expression of the minimum energy pathway (MEP) through such a barrier Φ(r,c). The MEP was therefore computed numerically from a steepest descent method. For c 0 =0.04, a process in three steps is exhibited (fig. 4). The cluster composition stays equal to the initial solute concentration up to R s =6. This means that nucleation only starts at R s = 6. For 6<R<8, the cluster composition sharply increases up to the expected value (c = 1). The change in size (R) related to this sharp enrichment of embryos is rather small. For R>8, the cluster grows without change of composition. The trajectory for c 0 =0.07, c 0 =0.03 and c 0 =0.02 are similar (Fig. 4). This non-classical model of nucleation gives rise to unanticipated effects. Surprisingly, nucleation is indeed observed to start for smaller radius R, when the supersaturation decreases (i.e. c 0 decreases). This is the consequence of the modification of the topology of the barrier Φ(r,c) (fig. 3) when changing the nominal composition. This modification, that is the result of computation, is complex and is almost impossible to foresee. This three-step nucleation process has already been found by Schmelzer et al. for η=0.7 [START_REF] Schmelzer | On the Dynamics of First-Order Phase Transitions in Multi-Component Systems[END_REF]. However our results show that this three-step process is not universal. Indeed, for smaller supersaturation (c 0 =0.01 and η=0.3 in Fig. 4), there is no more plateau before nucleation starts. The cluster composition increases sharply for c~0.46 and for R close to zero.

After this step, the composition slowly increases towards the equilibrium value.

We shall now examine the influence of the misfit on the minimum energy pathway for various values of the elastic energy parameter (k=0,1,2). For the same temperature, η=0.3 and nominal composition (c 0 =0.04, Fig. 5) the lattice misfit is shown to slow down the evolution of the nuclei composition with R. Increasing k increases the time required to reach the expected composition c=1 (k=0,1). However increasing further k (k=2, Fig. 5) reduces the onset size (R s ) for nucleation. A similar behaviour was found for c 0 =0.02 (Fig. 6). For increasing misfit, the size required for the embryo to reach its asymptotic composition increases.

Model versus experiments

There are several systems where experiments show that during the early stages of precipitation, nuclei do not have the expected composition. This is the case of Fe-Cu [47][48][49][50][51] and Fe-Cr [52] where temporal evolution of the nuclei composition has been observed. In boron implanted silicon where the expected precipitated phase has a different structure compared to the diamond lattice of silicon, very dilute boronenriched clusters (7-20 at.% of boron -1-5 nm in size), thought to be coherent with silicon, were evidenced [16,[53][54][55]. This suggests that the problem is general and not simply limited to the early stage of precipitation in metallic alloys. We shall apply the Diffusion mechanisms, particularly boron interactions with self-interstitials should be studied in more details in order to clarify the role of kinetics in the nucleation process in this system. Note that artefacts may also affect composition measured by atom probe tomography (APT). Local magnification effects may lead to underestimate the solute content in very small nuclei [57]. This can be particularly pronounced for SiB as boron is known to have a high evaporation field (defocusing effects in APT images [57]). Copper is indeed known to have a lower evaporation field compared to iron (focusing effects). Consequently, iron coming from the matrix may fall in the precipitate region in APT images. This may artificially increase its iron content (i.e. lower Cucontent). The influence of these effects has been studied theoretically and simulations predict them to be fairly small [57]. After thermal annealing or neutron or electronirradiation, a substantial fraction of Fe was detected in precipitates in the early stages, typically form 50 to 20 at.% (i.e. 50-80 at.% Cu for R=1-2 nm) [47][48][49][50][51]. Various interpretations accounting for the substantial amounts of Fe observed in Cu precipitates were proposed such as coherency stresses that might stabilize the precipitates or kinetic considerations [60]. We have applied the present model to the Fe-Cu system in order to establish whether thermodynamics considerations can account at least partly for the nucleation of clusters containing less than 100 at.% of Cu, the expected equilibrium value. We have considered the Fe-Cu system at 773K and for a nominal composition c 0 =0.015. The free energy curve, f(c), was obtained from Thermo-Calc® [START_REF]Thermo-Calc®[END_REF]. The The critical radius R* is thus close to 0.25 nm (=2σ/∆f with σ=0.28 J/m 2 derived from the nearest neighbour broken bond model [START_REF] Yang | [END_REF]). This time, this suggests that thermodynamically it is easier to form a critical nucleus with a much lower Cu-content than that of the equilibrium phase [47][48][49][50][51]. The minimum energy pathway during nucleation was computed using a steepest descent method along the nucleation barrier and is shown in Fig. 7. The black dashed line indicates the nucleation pathway in the present theory. This pathway passes through the saddle point (c=0.75, r=1). The nucleation barrier is overestimated in the classical nucleation theory (grey trajectory in Fig. 7 is a vertical line located at c = 1). It was found that the equilibrium composition (c ~ 1) is only reached when the radius approaches 1 nm (~1.5 R*). This indicates that even if nuclei are fairly dilute, they rapidly enrich up to the equilibrium. This behaviour is similar to the evolution of the nucleus composition simulated by Nagano et al. [63]. They found that during aging the nucleus composition sharply increases up to unity without any significant change in size. The particle growth then proceeds. This model predicts (for c 0 =0.015 and T=773 K) that for size larger than 1.5 R*, pure Cu-enriched nuclei should be observed. As the critical radius is small (< 0.5 nm), the present approach does not explain the large iron content observed in larger precipitates (typically R>2 nm) observed by atom probe methods in similar conditions [47][48][49][50][51].

Conclusion

This paper has dealt with the non-classical nucleation of coherent precipitates that have a solute content lower than the "equilibrium" phase. This issue not only concerns the early stages of coherent precipitation (e.g. FeCr) but also the nucleation of transient coherent phases preceding the precipitation of incoherent phases (e.g. FeCu, SiB). The generalized Gibbs model developed by Schmelzer et al. [35] has been Properties of critical nucleus (R*,c*) were derived from the saddle point in the R,c space for a regular solid solution (T/T c =0.3) and investigated as a function of the supersaturation. The model confirms that the classical nucleation theory is only valid for small supersaturation. Surprisingly, it was demonstrated that in the present theory the critical composition is independent on the interfacial coefficient that controls the dependence on the interfacial energy with the concentration gradient. The composition of critical nucleus was found to approach the matrix composition near the spinodal line.

The critical radius of nucleation was found to diverge to infinity not only near the solubility limit c α but also close to the spinodal line c s . The driving force and the nucleation barrier were both found to vanish near the spinodal line. The transition between the nucleation-growth and the spinodal regimes is found to be smooth in contrast to the classical nucleation theory.

The study of the influence of the elastic energy has shown non-anticipated effects. Paradoxically, increasing the elastic energy was found to increase the solute content in the critical nucleus. In spite of this non-classical behaviour, as expected, increasing the elastic energy increases both the nucleation barrier height and the critical radius. The minimum energy pathway was computed so that to predict the kinetic pathway of nucleation. We found that nucleation often proceeds with three step process.

Nuclei grow without change in composition and then they rapidly enrich until the

  at the interface. The nucleation barrier, proportional to σ 3 ,

  , the elastic energy related to misfit stress between phases was not considered despite of its crucial importance for precipitation in solids.

  Let us consider the nucleation of B-enriched precipitates in A-B binary solid solution (α). The contribution of the misfit elastic energy has been considered.

7 with

 7 γ% =E.δ 2 /(1-υ), with E the Young modulus, υ the Poisson coefficient and δ the misfit between pure phases (δ = δa/a, with a the average lattice parameter and δa = a Ba A ).

  coordinance and ε the energy of ordering classically expressed as a function of effective pair interaction energies between first nearest neighbour atoms (AA, BB, AB).

  β (right part of the binodal curve).

  2.b) whereas it reaches a maximum in the classical case. This result (Γ=0) was expected as c* tends toward c 0 when approaching the spinodal line. The work, W, required to form a critical nucleus (Fig.2c, with 16πσ 3 /3(ρk B T) 2 = 2.33 and 2σ/ρk B T = 12.7) is obviously lower than the classical nucleation barrier. W→0 when the initial solute concentration approaches the spinodal line. As a result, the critical radius tends to binodal or the spinodal lines are approached (Fig2.d). This confirms that the classical nucleation theory is only valid for small supersaturation for which cluster properties (composition and critical radius) are close to the classical theory and that the nucleation regime does not apply anymore when approaching the spinodal line.

  for greater c* when increasing the coherency misfit. This could appear as paradoxical as larger values of c* increase the elastic energy and as a result the nucleation barrier (equation (9)). However, in the present non-classical theory, the effective driving force also depends on c* and may win in the competition.Moreover, as the addition of the elastic energy corresponds to an increasing of T/T c , the minimal composition of the critical cluster, c m , increases. In contrast to the classical theory, the elastic energy influences the composition of embryo that in turn has a direct influence on the driving force. The chemical driving force Γ(c*) increases with the elastic energy because c* increases (Fig.

Fig. 3 (

 3 Fig.3(no elastic energy). The kinetic pathway has to pass through the saddle point (r=1

  present model to theses systems in order to highlight whether thermodynamic considerations (i.e. minimization of the nucleation barrier) can explain the nucleation of non equilibrium transient coherent clusters during the early stages of precipitation.Dilute boron-enriched clusters were observed in implanted silicon. Due to the quite different structure of the SiB 3 equilibrium phase compared to the diamond structure of the Si parent phase, the direct nucleation of the equilibrium phase appears difficult if not impossible. Instead, precipitation of this phase should proceeds through the formation of coherent isostructural boron-enriched clusters having the structure of the parent silicon phase. Experiments indeed showed the formation of dilute boron-typically 7-20 at.% B, 1-5 nm in size)[16]. The nucleation of such boron-enriched clusters at supersaturations and temperatures corresponding to the experimental conditions was modelled. The solubility limit of boron in silicon c α is close to 10 -4 at 873 K. In a regular solution model, the critical temperature is given by equation (12), T c ~4000K. This corresponds to T/T c ~0.22 (T=873 K). For a nominal composition c 0 of 1%, the minimum value c m of the critical composition is ~0.3. The computed value of the composition c* of the critical cluster is ~0.63 (driving force for nucleation (∆f) close to 7.5*10 8 J/m 3 ). The critical radius R* was found close to 0.5 nm (=2σ/∆f with σ=0.2 J/m 2[56]), a reasonable value compared to experiments. If now we take a rather high elastic energy (misfit of 10%[16]), critical clusters become almost pure in boron (c*~0.96). The effective driving force (i.e. taking into account elastic energy) was found lower (~6.3*10 8 J/m 3 ). The nucleation barrier therefore increases in comparison to the zero misfit case (as in Fig.2) but it remains lower than that computed for lower solute content with a misfit of 10%. Thermodynamics therefore cannot thoroughly explain the formation of such dilute nuclei (less than 10% for R~1-2 nm) observed during the early stages of precipitation in boron implanted silicon. Kinetics arguments should also be taken into account to explain the formation of dilute nuclei.

  The early stages of precipitation in binary Fe-Cr alloy for a nominal composition (c 0 ) below the spinodal line were recently investigated[52]. Using atom probe tomography, Novy et al. showed that the composition of the Cr-enriched precipitates (~0.62 for R~1-1.5 nm[52]) does not reach equilibrium (c β ~0.84) during the early stages of precipitation. This was modelled in the same conditions as this experimental work (i.e. T=773K, c 0 =0.2). Instead of a regular solution model, we have used the rather complex free energy function given byBonny et al. [see appendix of ref 58] so that to compute the driving force for nucleation as a function of the nucleus composition. This gives a solubility limit c α close to 15% at 773 K (in good agreement with experimental phase diagram) and the spinodal is c s ~ 34%. The critical nucleus composition c* was found equal to 0.78 (eq.8 neglecting the elastic energy). The addition of a low elastic energy (typically a misfit of 1%) does not change c* because it is negligible in comparison to the driving force. The critical nucleus composition c* get closer to the equilibrium state (c β ~0.84) than that observed by Novy et al. (c~0.62 for R~1-1.5 nm [52]). The corresponding driving force for nucleation (∆f) is close to 3*10 8 J/m 3 . The critical radius R* is thus close to 1 nm (=2σ/∆f with σ=0.15 J/m 2 [59]). Our model based only on thermodynamics considerations shows the good trend but again it does not fully explain such the Cr content in the early stages of precipitation in Fe-Cr. Again, kinetics effects related to diffusion should play a role.Precipitation of Cu (FCC) in the low temperature BCC phase of Fe is known to proceed via the nucleation of transient coherent embryo that have the same BCC structure as the parent Fe phase. The composition of such bcc Cu-Fe precipitates is still a controversial issue[47]. One source of this is that composition data estimated using APT may be biased because of the occurrence of already mentioned local magnification

corresponding solubility limit is close to 10 - 4 .

 104 The minimum value c m of the critical composition is ~0.4. Equation (8) (no misfit) gives a critical composition in copper, c*, equal to 0.75. As in Fe-Cr, the elastic energy addition (misfit of 1%) does not change drastically c*. The driving force for nucleation (∆f) is close to 2.2*10 9 J/m 3 (no misfit).

  to include the influence of the elastic energy originating from the lattice misfit between nuclei and the parent phase. Isotropic elastic energy related to misfit has been considered. The composition of critical nuclei was determined through the minimization of the nucleation barrier in the radius (R) and composition (c) space.

Fig. 1 Fig. 2

 12 Fig.1Schematic representation of the free energy curve for a system exhibiting isostructural decomposition.

Fig. 3 Fig. 4

 34 Fig.3The reduced nucleation barrier Φ as a function of the nucleus size (reduced radius r = R/R*) and composition (c) for η=0.3, and a nominal composition c 0 =0.04 (k=0). The saddle point is localized at r=1 and c= c*=0.456.

Fig. 5 Fig. 6 Fig. 7

 567 Fig.5Evolution of the solute content in nuclei (c) as a function of their size (R) for various elastic constant (k= 0, 1, 2 and η=0.3). The nominal composition was c 0 =0.04.

  Nagano et al. have simulated the composition of Cu particles formed during isothermal aging in bcc Fe-Cu alloys using the Cahn-Hilliard nucleation theory and non-linear diffusion equation [63]. They found that whereas Cu nuclei are almost pure Cu in low Cu content alloys (<1 at.% Cu), the Cu concentration in the nucleus is significantly less than unity in higher Cu content alloys (>1 at.% Cu). Results derived from the present theory are in good agreement with the work of Nagano et al. as we found c*=0.75 for a high Cu content alloys (1.5 at.%).

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgements The authors acknowledge C. Pareige, B. Radiguet (GPM, Rouen) and G. Martin (CEA, Saclay) for fruitful discussion. The first author also acknowledges N. Masquelier for his help in the use of the Thermo-Calc® software.


E Sarrouy

O Dessombz
email: olivier.dessombz@ec-lyon.fr

J J Sinou

Piecewise polynomial chaos expansion with an application to brake squeal of a linear brake system

Keywords: brake squeal noise, stochastic, friction, eigenvalue problem

This paper proposes numerical developments based on polynomial chaos (PC) expansions to process stochastic eigenvalue problems efficiently. These developments are applied to the problem of linear stability calculations for a simplified brake system: the stability of a finite element model of a brake is investigated when its friction coefficient or the contact stiffness are modeled as random parameters. Getting rid of the statistical point of view of the PC method but keeping the principle of a polynomial decomposition of eigenvalues and eigenvectors, the stochastic space is decomposed into several elements to realize a low degree piecewise polynomial approximation of these quantities. An approach relying on continuation principles is compared to the classical dichotomy method to build the partition. Moreover, a criterion for testing accuracy of the decomposition over each cell of the partition without requiring evaluation of exact eigenmodes is proposed and implemented.

Several random distributions are tested, including a uniform-like law for description of friction coefficient variation. Results are compared to Monte Carlo simulations so as to determine the method accuracy and efficiency. Some general rules relative to the influence of the friction coefficient or the contact stiffness are also inferred from these calculations.

Introduction

Brake squeal noise is a real problem for the automotive [START_REF] Ouyang | A methodology for the determination of dynamic instabilities in a car disc brake[END_REF][START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF] and railway [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] industry: it is a discomfort for passengers, generates high warranty costs for car manufacturers and causes problems to train operators that have to face claims of people living next to train stations. Indeed, this complex phenomenon occurs at low speeds. It is acknowledged to come from the vibration of the brake components with high and audible frequencies and high intensity (up to 130 dB). Many deterministic studies are available in the literature using different mechanisms to explain brake squeal and going from the simplest models with a few degrees of freedom only [START_REF] Hoffmann | A minimal model for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF] to large finite element models taking many components into account [START_REF] Cao | Linear eigenvalue analysis of the disc-brake squeal problem[END_REF]. Among the mechanisms put forward to explain this noise generation, one finds the decreasing friction coefficient with sliding speed, sprag-slip [START_REF] Spurr | A theory of brake squeal[END_REF], stick-slip [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal, and chaos-part i: Mechanics of contact and friction[END_REF][START_REF] Ibrahim | Friction-induced vibration, chatter, squeal, and chaos-part ii: Dynamics and modeling[END_REF] and self-excited vibration coming from mode coupling with a constant friction coefficient [START_REF] North | Disc brake squeal, a theoretical model[END_REF][START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF]. Besides the above mentioned references, an extensive review is available in [START_REF] Kinkaid | Automotive disc brake squeal[END_REF][START_REF] Ouyang | Numerical analysis of automotive disc brake squeal: a review[END_REF]. In this paper, a simple Coulomb friction law with a constant friction coefficient is used. Instability then occurs when two modes couple. This can be detected by analyzing the eigenvalues of the linear problem or the tangent one if non-linear effects are considered. In the latter case it is also possible to determine the limit cycles using frequency domain methods [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF][START_REF] Chevillot | Effects of damping on the speed of increase and amplitude of limit cycle for an aircraft braking system subjected to mode-coupling instability[END_REF] or to conduct a transient analysis [START_REF] Chern | Nonlinear brake squeal analysis[END_REF][START_REF] Sinou | Transient non-linear dynamic analysis of automotive disc brake squeal -on the need to consider both stability and non-linear analysis[END_REF]. In this work, we focus our attention on the stability analysis of a linear brake system with uncertainties. The stability of the equilibrium of the linear system is investigated via the Complex Eigenvalue Analysis (CEA) method. The eigenvalue problem has then to be solved. The objective of the present study is to propose a methodology and numerical developments based on polynomial chaos expansions (PCE) to process stochastic eigenvalue problems accurately and efficiently and to save time compared to Monte Carlo simulations. As no non-linear phenomena are considered, the prediction of squeal frequencies are limited. Nonetheless, it is a first step towards faster calculations for processing brake systems including uncertainties.

The present model is a linear finite element model (FEM) whose equilibrium stability is investigated through the corresponding eigenvalue problem. Instead of refining the model to get precise information on its vibratory behavior, some parameters are considered to be random. As pointed out by several papers [START_REF] Bergman | Influence of disc topography on generation of brake squeal[END_REF][START_REF] Oberst | Statistical analysis of brake squeal noise[END_REF], some parameters such as the friction coefficient are hard to characterize accurately for each brake and each braking action. It then seems appropriate to conduct a stochastic study to take into account this variability.

The main problems encountered for a stochastic study is the description of the random inputs and the numerical cost [START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF]. Indeed, a well known and robust method is the Monte Carlo simulation (MCS). It consists in the evaluation of quantities of interest (for example here, the eigenvalues) for a large sample of uncertain inputs values. Statistical estimators (such as the mean, standard deviation, etc.) can then be rebuilt (see for example [START_REF] Culla | Uncertainty model for contact instability prediction[END_REF] for MCS applied to a brake model). For systems with more than a few degrees of freedom, this operation can have a prohibitive cost in terms of computation time even though recent developments relying in particular on parallel processing [START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF] let MCS remain an attractive method. Apart from MCS, three other classes of methods were developed [START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF]. First the perturbation methods are based on Taylor series expansion on a set of zero mean random variables [START_REF] Liu | Random field finite elements[END_REF][START_REF] Hien | Finite element analysis based on stochastic hamilton variational principle[END_REF]. The series are usually truncated to first or second order to avoid difficulties linked to the evaluation of higher order derivatives. Hence they are mostly used in the case of small variation of random inputs. It is to be noted that recent developments allow the computation of complex eigenvalues [START_REF] Adhikari | Complex modes in stochastic systems[END_REF][START_REF] Rahman | Stochastic dynamic systems with complex-valued eigensolutions[END_REF]. The second family of methods relies on the so-called Polynomial Chaos Expansions and was first introduced in the field of structural dynamics by Ghanem and Spanos [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. It basically consists in the decomposition of quantities of interest on an orthogonal basis of polynomials (see for example [START_REF] Sepahvand | Uncertainty quantification in stochastic systems using polynomial chaos expansion[END_REF] for a detailed presentation of such methods). The procedure proposed in this paper derives from PCE and its further developments. A new method to build the stochastic space partition is proposed, based on continuation method principles. Moreover, an efficient accuracy criterion dedicated to the eigenvalue problem is developed to test the current polynomial decomposition over each element of the partition without having to compare it to any exact result. These several improvements aim at providing a method that returns accurate results in a lesser time. Section 2 exposes the method in detail. Finally, Soize [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Murthy | Nonparametric stochastic modeling of uncertainty in rotordynamics-part i: Formulation[END_REF] developed a non-parametric approach which does not require identification of the uncertain local parameters. A deeper overview of these several methods is available in successive review papers by Ibrahim et al. [START_REF] Ibrahim | Structural dynamics with parameter uncertainties[END_REF][START_REF] Manohar | Progress in structural dynamics with stochastic parameter variations: 1987-1998[END_REF] and Schüeller et al. [START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF][START_REF] Schuëller | A state-of-the-art report on computational stochastic mechanics[END_REF][START_REF] Schuëller | Computational stochastic mechanics -recent advances[END_REF].

This paper is organized as follows: first, Section 2 exposes the theoretical background. Then the stochastic study of a FEM brake model is performed using the proposed method for several random laws. Accuracy and efficiency of the method compared to MCS is evaluated and the stochastic results about system stability are commented.

Theoretical background used to process uncertainty

This section is devoted to the presentation of the mathematical background of the method proposed in order to handle the eigenvalue problem of uncertain systems. First, quick recalls are proposed about the stability study of equilibria through the tangent eigenvalue problem. Then, the general idea of using an approximation of quantities of interest by polynomial decomposition is studied. Finally, a refinement of this approach that enables the study of uncertain quantities with strong variations is presented. The last part is dedicated to the algorithm that organizes theses different points.

Deterministic study of stability through the eigenvalue problem

The deterministic equation for a general linear structure is:

M q + D q + Kq = f e (1)
where q, q and q denote the vector of degrees of freedom (dofs) and its first and second derivatives with respect to time. M, D and K are the mass, damping and stiffness matrices respectively and f e is the vector of excitation forces. The static equilibrium for such a formulation is solution of Kq s = f e . Its stability depends on real parts of eigenvalues of the following problem [START_REF] Géradin | Théorie des Vibrations, Application à la dynamique des structures[END_REF]:

Au k = λ k Bu k , with A = 0 K -K -D and B = K 0 0 M (2)
Computing the 2n eigenvalues λ k and the matching eigenvectors u k solutions of Eq. ( 2) (where n denotes the system size), one can determine the stability of the equilibrium: if there exists at least one eigenvalue with strictly positive real part, then it is unstable, otherwise, this is a stable equilibrium position. Moreover, in order to define the eigenvectors u k uniquely, one will use normalization with respect to B matrix:

u T k Bu k = 1 (3)

Introducing uncertainties

If some parts of the structure have uncertain properties, they can be rendered by adding uncertain terms to the dynamic system (1):

M + M q + D + D q + K + K q = f e (4)
Random contributions are denoted with tilde notation. It is assumed that they have zero mean and depend on a set of random parameters denoted ξ.

The stability study of the static equilibria of Eq. ( 4) is then achieved through the following eigenvalue problem

A + Ã ũk = λk B + B ũk (5)
with normalization equation

ũT k B + B ũk = 1 (6)
where ũk and λk denote the k-th stochastic eigenvector and eigenvalue respectively.

The proposed method to compute these stochastic eigenmodes uses a decomposition of uncertain quantities on a polynomial basis. Depending on the law followed by the set of random parameters ξ, some basis are optimal considering the convergence rate with polynomial degree [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF]. For example Hermite polynomials are most adapted to normal law while Legendre polynomials suit best the uniform law. One of the most popular methods based on such a polynomial decomposition is the Polynomial Chaos expansion which uses a basis of Hermite multivariate polynomials [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Ghosh | Analysis of eigenvalues and modal interaction of stochastic systems[END_REF]. A previous work by the authors addresses the dispersion of a finite element rotor modes using such a basis [START_REF] Sarrouy | Stochastic analysis of the eigenvalue problem for mechanical systems using polynomial chaos expansion -application to a finite element rotor[END_REF]. Recently, this method has been extended to the polynomials of the Askey scheme leading to the so-called "generalized Polynomial Chaos" [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF]. Note that the expression "polynomial chaos" refers to the theory following Wiener's developments [START_REF] Wiener | The homogeneous chaos[END_REF] and is to be distinguished from the "deterministic chaos" [START_REF] Oberst | Chaos in brake squeal noise[END_REF]. The rest of this section is dedicated to the presentation of the method adapted to the stochastic normalized eigenvalue problem [START_REF] Cao | Linear eigenvalue analysis of the disc-brake squeal problem[END_REF] and [START_REF] Spurr | A theory of brake squeal[END_REF].

To this end, let us denote ψ ψ n , 1 ≤ n ≤ N the orthogonal polynomial basis with N elements used for uncertain quantities decomposition and < ψ ψ i , ψ ψ j > the associated scalar product. First, matrices introducing uncertainty are decomposed as follows, using a Karhunen-Loève decomposition or assuming them so regarding experimental results [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF]:

X = N X n=2 X n ψ ψ n (ξ) (7)
where X can denote M, D or K. ψ ψ 1 will conventionally denote the constant polynomial equal to 1 (associated with deterministic component); this implies that polynomials with index greater or equal to 2 are zero mean. Then stochastic eigenvalues and eigenvectors are decomposed on both the deterministic eigenmodes and the polynomial basis:

λk = λ k N s=1 (k) a s + j (k) b s ψ ψ s (ξ) (8) ũk = p∈P k (k) γp (ξ) + j (k) μp (ξ) u p (9) with (k) γp (ξ)
= N n=1 (k) γ n p ψ ψ n (ξ) and (k) μp (ξ) = N n=1 (k) µ n p ψ ψ n (ξ) (10)
giving ũk = p∈P k        N n=1 (k) γ n p + j (k) µ n p ψ ψ n (ξ)        u p ( 11
)
where the coefficients (k) γ n p , (k) µ n p , (k) a s and (k) b s are real and j denotes the imaginary unit (j 2 = -1). P k is the set of deterministic modes used for the k-th stochastic mode decomposition. It contains P k elements. This decomposition is a generalization of what is proposed in [START_REF] Dessombz | Analyse dynamique de structures comportant des paramètres incertains (Dynamic analysis of structures with uncertain parameters[END_REF]. The use of complex weights (k) γ n p + j (k) µ n p and (k) a s + j (k) b s lets real and imaginary parts of stochastic eigenvectors and eigenvalues evolve independently. The selection of a few deterministic modes included in P k instead of using all the deterministic modes (P k = 2n) can lead to a great computation time reduction and is addressed in Section 3.3. To get the final system of equations that leads to the unknowns (k) γ n p , (k) µ n p , (k) a s and (k) b s , equations ( 5) and ( 6) are first projected onto the polynomial basis ψ ψ n , 1 ≤ n ≤ N using the dedicated scalar product. This gives the set of equations ( 12) and [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF]

for 1 ≤ i ≤ N. p∈P k N A n A =1 (A n A u p )        N n=1 ( (k) γ n p + j (k) µ n p ) < ψ ψ n A ψ ψ n , ψ ψ i >        -λ k p∈P k N B n B =1 (B n B u p )×         N n,s=1 ( (k) γ n p + j (k) µ n p )( (k) a s + j (k) b s ) < ψ ψ n B ψ ψ s ψ ψ n , ψ ψ i >         = 0 (12)
and

p 1 ,p 2 ∈P k N B n B =1 (u T p 1 B n B u p 2 )×         N n 1 ,n 2 =1 ( (k) γ n 1 p 1 + j (k) µ n 1 p 1 )( (k) γ n 2 p 2 + j (k) µ n 2 p 2 ) < ψ ψ n B ψ ψ n 1 ψ ψ n 2 , ψ ψ i >         -1 < ψ ψ 1 , ψ ψ i >= 0 ( 13
)
where A n A and B n B are the matrices used to decompose random matrices à and B respectively, in accordance with Eq. ( 7) notation.

-1 -0. Separating real and imaginary parts of these equations, one gets 2N(2n+1) real scalar equations for 2N(P k +1) real unknowns. In case of a deterministic mode set P k which does not include all the deterministic modes (P k 2n), this system is overconstrained. To get a square system, the final step is to pre-multiply Eqs. [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal: a review[END_REF] by the transpose of the matrix containing the deterministic modes in P k . This finally leads to a square non-linear (quadratic) system with size 2N(P k + 1) that can be processed using a dedicated non-linear solver.

This procedure lets us approximate eigenmodes λk and ũk by a complex sum of multivariate polynomials ψ ψ n . The usual statistical estimators such as mean and standard deviation can then be quickly evaluated over a large sample of ξ points by using formulas [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal, and chaos-part ii: Dynamics and modeling[END_REF] and [START_REF] Kinkaid | Automotive disc brake squeal[END_REF]. This is equivalent to a Monte Carlo simulation but without the resolution of an eigenvalue problem for each realization and is much faster. As this is a polynomial approximation, one needs to use high degrees to approximate quantities with strong variations accurately. For example, if an equilibrium becomes unstable over the input range under study, the first derivative of some eigenvalue may not be continuous; this makes the polynomial approximation very poorly appropriate to describe the eigenvalue evolution. From a more general point of view, studying changes in stability, it is obvious that different evolutions will be observed over the different stable and unstable ranges. This is illustrated on Fig. 1. On this figure, the target curve has a discontinuous derivative and y axis matches for example the evolution of the imaginary part of an eigenvalue whose mode undergoes a change in stability for ξ = -0.2. The several polynomial approximations using different polynomial degrees displayed in Fig. 1 show that a high degree (pane (c) degree 10) is required to describe the change in slope accurately. Decompositions using a polynomial basis with such a high degree generate very large systems ( 12) and ( 13) to solve. To overcome this problem, one may use a piecewise polynomial approximation. The next section is devoted to the introduction of this piecewise approach into the general method.

Working with a partitioned stochastic space

The idea of using piecewise polynomial approximation of uncertain quantities has already been developed in the field of fluid mechanics [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Le Meitour | Prediction of stochastic limit cycle oscillations using an adaptive polynomial chaos method[END_REF] for the study of limit cycles or equilibria. The subsequent method is referred to as Multi-Element generalized Polynomial Chaos (MEgPC). These studies use a decomposition of the random space to which the random variable vector ξ belongs. A polynomial approximation is then computed over each element of this decomposition. It is generally developed considering uniform laws for the random variable ξ which leads to using Legendre polynomials for the approximation. This restriction to a uniform law can be alleviated as non-uniform random variables can be approximated by a sum of Legendre polynomials of uniform random variables [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF].

In this section, the decomposition of the random space and its consequences on the stochastic eigenmodes approximation are exposed with a non-probabilistic point of view: it is presented as a piecewise polynomial fitting of eigenvalues and eigenvectors components depending on a vector of Q real variables grouped in vector ξ. The aim of the procedure is then to get the expressions of eigenmodes using a decomposition on a Legendre orthogonal polynomial basis on each element of the partition. Different statistical laws for uncertain input description are then handled by adapting the sample used for the Monte Carlo simulation. This MCS is expected to be realized quickly as the eigenmodes are evaluated using the piecewise polynomial approximation rather than by using a direct resolution of the eigenvalue problem. One could object that the Legendre polynomial basis is not the one that suits best the input distribution (for example Gamma distribution) and hence, that the weight function used in the associated scalar product does not put stress on the most probable values. This is true except if the resulting description is accurate enough whatever the chosen point may be.

First the description of a quantity depending on ξ belonging to S i , a bounded interval of IR Q using a basis of Legendre multivariate polynomials depending on a variable ζ belonging to I Q is exposed. Then, the problem of switching from such a description over S i to the description over S j , included in S i is addressed. The global scheme using these fundamental steps is presented in Section 2.4.

Let ψ ψ n (ζ), 1 ≤ n ≤ N be the Legendre multivariate polynomials basis with degree less or equal to D in the case of Q variables, defined over

I Q , I = [α, β]. The scalar product associated is < f, g >= I Q f (ζ)g(ζ)w(ζ)dζ ( 14
)
where w = 1 denotes the weight function associated to this polynomials family.

Let S i be a bounded interval of IR Q and denote ν i the diffeomorphism mapping S i into I Q :

ν i : S i → I Q ξ → ν i (ξ) = ζ (15)
If

S i = [a (i) 1 , b (i) 1 ] × • • • × [a (i) Q , b (i) Q ], one will use ∀k ∈ [[1, Q]], ζ k = αb (i) k -βa (i) k b (i) k -a (i) k + β -α b (i) k -a (i) k ξ k as ν i definition.
If f is a scalar function of ξ defined over S and {S i } i=1...N S is a partition of S , its value over S i will then be computed using:

f (ξ) = N n=1 f (i) n ψ ψ n (ν i (ξ)) (16)
The f (i) n are the coefficients of the decomposition of f on the polynomial basis ψ ψ n over S i . Considering the eigenvalue problem at stake, on will solve the following problem on each subset S i :

Ã(i) ũ(i) k = λ(i) k B(i) ũ(i) k (17) ũ(i) k T B(i) ũ(i) k = 1 ( 18
)
with the following notations for each quantity involved:

Ã(i) (ξ) = N A n=1 A (i) n ψ ψ n (ν i (ξ)) (19)
B(i) (ξ) = N B n=1 B (i) n ψ ψ n (ν i (ξ)) (20)
λ(i) k (ξ) = λ k N n=1 (k) a (i) n + j (k) b (i) n ψ ψ n (ν i (ξ)) (21)
ũ(i) k (ξ) = p∈P k        N n=1 (k) λ p(i) n + j (k) µ p(i) n ψ ψ n (ν i (ξ))        u p (22)
To get the coefficients A (i) n and B (i) n , one must be able to deduce the decomposition of any quantity over an interval S j included in S i knowing the decomposition over S i . This procedure is explained in the case of a scalar quantity f . It is easily adaptable to non-scalar quantities. Let S j be an interval included into S i and g the restriction of f over this interval. g decomposition over the ψ ψ n family is denoted:

g(ξ) = N n=1 f ( j) n ψ ψ n (ν j (ξ)) (23)
Knowing the coefficients f (i) n , one can compute the f ( j) n coefficients in the following manner: first, consider ξ ∈ S j and look for decomposition coefficients χ n,k such that

ψ ψ n (ν i (ξ)) = N k=1 χ n,k ψ ψ k (ν j (ξ)) (24)
The χ n,k coefficients are found by solving for each 1 ≤ n ≤ N the system:

            ψ ψ 1 (ν j (ξ 1 )). . . ψ ψ N (ν j (ξ 1 )) . . . ψ ψ 1 (ν j (ξ N )). . . ψ ψ N (ν j (ξ N ))                          χ n,1 . . . χ n,N              =              ψ ψ n (ν i (ξ 1 )) . . . ψ ψ n (ν i (ξ N ))              , ξ p ∈ S j ( 25
)
That is, N linear systems with size N are solved to establish the χ matrix such that:

∀ξ ∈ S j ,              ψ ψ 1 (ν i (ξ)) . . . ψ ψ N (ν i (ξ))              = χ              ψ ψ 1 (ν j (ξ)) . . . ψ ψ N (ν j (ξ))              (26)
Finally, the f ( j) n coefficients are obtained from the f (i) n using

f ( j) 1 , . . . , f ( j) N = f (i) 1 , . . . , f (i) N χ (27)

The global algorithm: two possible strategies

To present the global algorithm, we assume that the random inputs are described using Q uniform random variables ξ 1 , . . . , ξ Q that vary over S , a bounded interval of IR Q . The aim of the algorithm is to return a collection of polynomial approximations of eigenmodes using a Legendre basis with N elements over a partition of S whose elements are denoted S i ; elements of the partition will be referred to as "cells". Over each S i , one then has to determine the coefficients (k)

λ p(i) n , (k) µ p(i) n , (k)
a (i) s and (k) b (i) s .
For notation convenience, these 2N(P k + 1) scalar coefficients are gathered in a vector denoted x (i) . Two strategies to process the stochastic space are deployed and will be compared in Section 3: the first one relies on the dichotomy method and is close to the method proposed in [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Le Meitour | Prediction of stochastic limit cycle oscillations using an adaptive polynomial chaos method[END_REF]. This method is exposed in Section 2.4.1. The second one invokes continuation principles and is presented in Section 2.4.2.

That said, the two strategies share a common scheme presented in Fig. 2. First, the algorithm is initialized with the set of matrices A n and B n that define the problem over S and a first guess for decomposition coefficients x 0 ; this vector can for example be chosen such that the processed stochastic eigenmode matches the deterministic mode over the whole S interval. Then, the general loop is applied: first, if all the space S has been processed, algorithm stops. Otherwise, a cell S i is defined in the remaining space as well as a startpoint x (i) i , depending on the chosen method. The collection of matrices A (i) n and B (i) n defining the problem over the current cell are evaluated using ( 26) and ( 27) and the problem ( 12) and ( 13) is then solved, returning a new set of coefficients x (i) f . Accuracy of the decomposition is tested as detailed in Section 2.4.3. If it satisfies the criterion or the criterion is not satisfied but no more cell can be created (a maximum computation time, a maximum number of cells or a minimum cell width has been reached), each element is stored (that is the current cell definition, eigenmode decomposition coefficients, matrices decomposition and accuracy criterion result). Otherwise, the current cell S i is refined accordingly to the chosen method.

First method: Dichotomy

The simplest way to refine a cell over which the decomposition is not satisfying is to split it in half in one direction to generate two smaller sub-cells. Then, a startpoint has to be chosen for future processing of each of these sub-cells; different choices are possible for this startpoint. First, one can always use the same startpoint x 0 that is provided by the user to initialize the algorithm. Another possibility is to convert the final decomposition x (i) f over S i into equivalent decompositions over each of the sub-cell using formulas [START_REF] Sepahvand | Uncertainty quantification in stochastic systems using polynomial chaos expansion[END_REF] and [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]. Finally, one can extrapolate a decomposition from neighboring cells: considering a one dimensional problem, once the decomposition over one cell is validated, the startpoint for each of the two neighbors is the decomposition that matches constant values for eigenmodes over each sub-cell. The value is equal to the value of eigenmodes at the validated cell border. Fig. 3 illustrates this procedure: once the decomposition over S i is validated, each x (i) component x j of neighboring cells S left ans S right is initialized with the value at corresponding border of S i . This procedure giving better results than the two previous ones mentioned above, it is the one that is used for numerical applications of Section 3.

Second method: Continuation

Another point of view is to consider that the sought values are continuous from one cell to another and hence, to apply some continuations principles (see [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF] for the continuation methods). Hence, in one dimension, starting from the center of the initial cell, new cells are built and processed by having ξ varying continuously to the right and then from the center to the left. Once a cell has a proper decomposition (satisfying the accuracy criterion), another cell -contiguous -is computed: its width depends on the current cell width and various parameters such as the number of iterations used to converge to the decomposition over this cell. If only a few iterations were used, the new cell width is increased, if the number of iterations exceeds a certain threshold, then the width of the next cell is decreased; otherwise, it is unchanged. When the decomposition over the current cell is not satisfying, the cell width is decreased too. The startpoint over each cell processed is computed using the contiguous cell decomposition which is extrapolated over this new cell. This method, compared to the one based on dichotomy injects more information at each cell by using extrapolation of the contiguous decomposition and using a more appropriate size for each cell. Its efficiency is tested in Section 3.

Testing the accuracy of the decomposition

Once the solution x (i) f has been evaluated for cell S i , one needs to decide whether the polynomial approximation is satisfying or not. That is, one wants to test if it is accurate enough over S i . As comparing value of eigenmodes returned by the approximation and a direct computation for several points over the interval would cost too much, one proposes an original criterion based on the Rayleigh ratio that is very cheap to evaluate.

Combining Eqs. ( 17) and ( 18), one gets the Rayleigh ratio

∀ξ ∈ S i , ũ(i) k T Ã(i) ũ(i) k ũ(i) k T B(i) ũ(i) k = λ(i) k ( 28
)
As the equations are satisfied in a weak sense only (only their projection over the polynomial basis are considered), this ratio may not be verified exactly. This provides a criterion to estimate the quality of the current decomposition by testing

ǫ k (ξ) =           ũ(i) k T Ã(i) ũ(i) k ũ(i) k T B(i) ũ(i) k           / λ(i) k -1 (29)
for several points ξ in S i .

Application to a finite element model for brake squeal

This section is devoted to the application of the proposed method to a brake model in order to conduct a stochastic study of brake squeal noise. First, the considered model is exposed along with the uncertain configurations. Then, the accuracy of the method using successively a dichotomy and a continuation approach is evaluated through the comparison of their results to direct resolutions. Finally, histograms are built, based on these results in order to support the stochastic study of the system and demonstrate its efficiency compared to Monte Carlo simulations.

Brake model and uncertain configurations

The model used is derived from the finite element model described in [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF][START_REF] Sinou | Study of the nonlinear stationary dynamic of single and multi-instabilities for disk brake squeal[END_REF]. It consists of a rotating disc and a fixed pad as displayed in Fig. 4. The disc and the pad are modeled using the finite element method. Each model is reduced using a Craig and Bampton technique: 4 contact nodes are kept on each part (pad and disc) as well as 4 dofs (degrees of freedom) on the backplate of the pad for application of hydraulic pressure; the first 30 eigenmodes are retained on each part. Such a reduction may degrade the pertinence of the results regarding the real behavior of the brake system. This problem should be addressed carefully in the frame of a study which aims at getting quantitative results. The current model, used in previously cited work, was provided as is to the authors. The assembled model has then 88 dofs and the set of equations can be written

M q + D q + Kq + K c q = f e ( 30
)
where M is the mass matrix, K is the structural stiffness matrix, D is the damping matrix with a Rayleigh expression D = αM + βK, K c takes into account contact forces acting between the pad and the disc and f e stands for the force due to brake pressure applied on the pad. K c is a non-symmetric matrix that couples the two subsystems -the pad and the disc. It consists in the addition of two matrices: a first one, symmetric, which represents the normal force at the contact surface and which is proportional to the contact stiffness denoted k c . The second matrix that contributes to K c is non-symmetric and stands for the tangent force at the contact interface; it is proportional to the product of the friction coefficient -denoted µ -and k c . The contact follows a Coulomb friction law with constant friction coefficient µ and assumption of permanent sliding is made; more complex phenomenon such as stick-slip or sprag-slip are then put aside. Normal and tangential contact forces applied by the pad over the disc are opposite to forces applied by the disc on the pad and can be formulated as follows for contact node i, 1 ≤ i ≤ 4:

F i,normal = k c (U i,Pad -U i,Disc ) (31)
F i,tangent = µF i,normal sign(V i,rel ) (32)
where U i,Pad/Dic and V i,rel denotes the normal displacement of node i for the pad or the disc and the relative velocity between the pad and the disc at contact node i respectively. This coupling of the two substructures affects both lower and higher frequency modes as shown by Fig. 5. The upper part (a) of this figure displays the absolute difference between the Modal Assurance Criterion (MAC, [START_REF] Allemang | The modal assurance criterion-twenty years of use and abuse[END_REF]) of the undamped uncoupled system compared to itself and the MAC of the undamped coupled system compared to the uncoupled eigenvectors (see Eq. ( 33)).

∆MAC = MAC(V uncoupled , V coupled ) -MAC(V uncoupled , V uncoupled ) (33)
Higher values of this difference let us detect the modes that are most affected by the coupling of pad and disc structures. One considers here that the j-th mode is significantly modified if one of the terms of the j-th column of ∆MAC is greater than 0.5. The lower part (b) then outlines eigenfrequencies of the coupled system that match these modes using red dots. Values used for this comparison are k c = 1 • 10 7 N.m -1 and µ = 0.25. As recalled in Section 2.1, the stability of any static equilibrium of Eq. ( 30) is determined by the following tangent problem:

M r + D ṙ + (K + K c )r = 0 (34)
One will then focus on the eigenvalues of this tangent system while including some uncertain parameters. Two nominal configurations are studied. A first one with a contact stiffness k c = 1 • 10 7 N.m -1 and a second one with a contact stiffness twice larger (k c = 2 • 10 7 N.m -1 ). These two configurations correspond to a tangent linear contact state obtained by experimental tests of pad compression [START_REF] Sinou | Study of the nonlinear stationary dynamic of single and multi-instabilities for disk brake squeal[END_REF]. For each configuration, the variation of the friction coefficient µ or the contact stiffness k c is evaluated. Uncertainty is introduced as follows:

µ(ξ) = µ 0 (1 + ξ∆µ) (35)
with µ 0 = 0.25 and ∆µ = 0.5, and

k c (ξ) = k c0 (1 + ξ∆k c ) (36)
with k c0 = 1 • 10 7 N.m -1 or k c0 = 2 • 10 7 N.m -1 and ∆k c = 0.5. ξ is a random variable following a uniform law on [-1, 1]. This means that for each configuration a variation of ±50% of the parameter µ or k c is examined. These large variation ranges will let us test the robustness of the method. Unstable modes evolution along with the random variable ξ are displayed in Fig. 6 for configuration 2, case 2 (k c varies). For each mode, both real and imaginary parts are plotted. This figure shows that the evolution of the studied modes can be quite complex and non-regular over the studied range. Configuration 1 and case 1 (when µ varies) would produce similar figures. Table 1 lists unstable areas per mode in the case when µ varies for each configuration and case.

Algorithm accuracy

The two versions of the procedure exposed in Section 2.4 were applied to each configuration for each case. Unstable modes were approximated one by one, using a Legendre polynomial basis with maximum degree D = 2. As suggested in Section 2.3, only a few deterministic modes were retained to describe the stochastic modes. Based on the results of Fig. 5, only the most impacted modes were kept (modes 1 to 34, 38 to 41 and 65 to 88 as well as their complex conjugates). This leads to a set P k with P k = 124 modes out of 176 for each unstable stochastic mode ũk . Accuracy is evaluated by comparing results obtained using a direct evaluation and approximated ones over 201 points equally spaced over the range [-1,1] of the random variable ξ. To report the accuracy of the method, global indicators are proposed instead of local ones that would generate a long list of values considering the total number of unstable modes that have been approximated. One will focus on the relative error over the real and imaginary parts of eigenvalues, the error on the sign of the real part, the error on the Hopf bifurcation point localization (point at which the real part of an eigenvalue crosses the abscissa axis having a non-null imaginary part) and finally the error over the eigenvectors even if these are quantities of lesser interest for this study. The formulas for the first two relative error estimates are:

e Reλ = Re(λ PC ) -Re(λ MC ) Re(λ MC ) × 100 (37)
e Imλ = Im(λ PC ) -Im(λ MC ) Im(λ MC ) × 100 [START_REF] Wiener | The homogeneous chaos[END_REF] where PC and MC superscripts respectively denote results obtained using the approximation and a direct evaluation respectively. The error on the sign of the real part of eigenvalues is a binary indicator telling whether the sign is the same for both approximated and direct evaluation or not. The error on the Hopf bifurcation point localization is the distance between ξ points for which the bifurcation occurs. This distance is to be compared to the width of the interval which is 2. Finally, the relative error on eigenvectors is evaluated using the following formula which returns a scalar indicator for each ξ point considered and each eigenvector:

e φ = max |u PC | -|u MC | max |u MC | × 100 (39)
Rather than displaying these indicators over the [-1,1] range, one will give the percentage of points that are below 2% and 5% of error in case of relative errors on eigenvalues, below 5% for relative error on eigenvectors, the percentage of points for which an error on the sign of the real part is committed and finally the distance between ξ points at bifurcations. Finally, these results are summed up in Table 2 for each configuration and case by returning the worst result among the several modes. Obviously, both versions of the proposed method return very accurate results: stability is well evaluated (the percentage of points whose real part sign is wrong, e sign(Reλ) , is lower than 0.5) and Hopf bifurcation points are properly identified with errors e ξHopf lower than 1 • 10 -3 for an interval width equal to 2. Moreover, errors on the eigenfrequencies e Imλ are very limited too, all of the approximated values being equal to the exact ones more or less 5%. Real parts of eigenvalues are also well estimated, particularly the one returned by the algorithm relying on continuation principles. Finally, it is worth to be noted that eigenvectors are estimated with a satisfying accuracy by the second method (based on continuation principles) especially, given that all the tested points present an overall error smaller than 5%. If these quantities are not of the utmost importance in the context of a linear stability study, they may be useful when one aims at finding the limit cycle, providing a useful direction in which to look for it or, to analyze variations of mode shapes whether the system undergoes a bifurcation or not. Considering this first set of results, both versions of the proposed method prove their ability to estimate eigenvalues of a uncertain system accurately, with a slight advantage for the version invoking continuation principles with a variable cell size and an extrapolation of neighboring cell result used as a startpoint.

Stochastic study and comments -Method efficiency

Now that the method has been tested and showed a satisfying accuracy, this section is devoted to the stochastic analysis of the brake system. As one of the advantages of this general method is to provide a decomposition accurate whatever the input random law is, the decompositions of stochastic eigenmodes are used to test different random input laws for each case: based on experimental work [START_REF] Oberst | Statistical analysis of brake squeal noise[END_REF] a uniform-like input law was designed to model friction coefficient variation. The following probability density function (pdf) over [-1, 1] was chosen to tend to reproduce the shape of experimental data available in this work:

         f (ξ) = -10.0446ξ 2 -15.9375ξ -5.8929 if -1 ≤ ξ < -0.8 f (ξ) = 0.1339(ξ + 0.8) + 0.4286 if -0.8 ≤ ξ < 0.8 f (ξ) = -16.7411ξ 2 + 26.9196ξ -10.1786 if 0.8 ≤ ξ ≤ 1 (40)
The continuous pdf is depicted in Fig. 7 (a), superimposed to the histogram describing the sample used for the simulation. This law is almost uniform, with a greater probability for higher values to occur. The subsequent results are compared to the ones obtained via a truncated Gaussian input law whose sample histogram is provided in Fig. 8:

f (ξ) = 1 √ 2π × 0.3 2 e -ξ 2 2×0.3 2 , -1 ≤ ξ ≤ 1 (41)
For case 2, when the contact stiffness k c is uncertain, a uniform law, Fig. 9 (a.1) and (a.2), is compared to a truncated normal law [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF] depicted in Fig. 10 ( ): first, the percentage of stable and unstable occurrences with specification of the number of unstable modes is displayed. This lets us analyze the propensity to squeal and the difficulty to prevent squealing (related to the number of unstable modes to stabilize). Second, an histogram depicting the unstable frequencies (that is the imaginary part of eigenvalues with positive real part) is provided. Finally, a third histogram displays the repartition of real parts of unstable modes. This gives information about the instability growth rate. The results displayed in Figs. 7, 8, 9 and 10 come from the second version of the method, based on continuation principles which was proven to be more accurate. Nonetheless, the other version of the method was applied too so as to provide a comparison of numerical resources consumption which is addressed in the last part of this section.

Considering the first case relative to friction coefficient µ variation for both configurations and two different input laws whose results are depicted in Figs. 7 and8, one can first denote that results notably vary when switching from a uniform-like input law to a truncated Gaussian one: while the simulation based on a truncated Gaussian law parameter variation returns a lower number of simultaneously unstable modes, it also decreases the number of stable occurrences for both configurations. Moreover, unstable frequencies in the range [START_REF] Dessombz | Analyse dynamique de structures comportant des paramètres incertains (Dynamic analysis of structures with uncertain parameters[END_REF]50] kHz that appear in the case of simulations using a uniform-like law for the second configuration are not represented when considering a truncated normal distribution. Similar remarks can be inferred from comparison of the second case simulations (when k c varies) relying on the one hand on a uniform law (Fig. 9) and on the other hand on a truncated Gaussian law (Fig. 10): the truncated Gaussian law returns more pessimistic results in terms of stability with respectively 0.6 % and 22 % of stable occurrences for configurations 1 and 2 to be compared to 7 % and 30 % in the case when a uniform distribution is used. These comparisons of results based on the same decomposition of uncertain eigenvalues for different random laws emphasize the interest of the method that lets us reuse the decompositions computed once for all and hence save time compared to Monte Carlo simulations as will be discussed at the end of this section. Indeed, results returned when using several input laws are widely different and would merely not lead to the same optimization objectives in an engineering process. As the description of a random parameter such as friction coefficient may not be easy, this ability to quickly simulate a large panel of laws, including non-classical ones (such as the proposed uniform-like law) is a non-negligible advantage.

Comparing the results obtained for configuration 1 (k c = 1 • 10 7 N.m -1 ) and configuration 2 (k c = 2 • 10 7 N.m -1 ) for the same random law leads to the conclusion that a more stiff contact leads to a greater number of simultaneously unstable modes: up to 7 modes can become unstable simultaneously for configuration 2 (Figs. 7 (c) and 8 (c)) whereas a maximum of 4 modes is observed in the first configuration (Figs. 7 (b) and 8 (b)) when µ varies and these are 3 (Figs. 9 (c) and 10 (c)) against 2 (Figs. 9 (b) and 10 (b)) unstable modes for the second case studying variations of k c . The number of stable occurrences merely increases with contact stiffness turning, for the first processed case, from 12.5 % to 24 % in the case of a uniform-like law for µ variation and from 7.6 % to 9.8 % in the case of a truncated Gaussian law. It is to be noted that a greater stiffness generates also a wider spectral range with higher frequencies (see "Unstable frequencies" histograms from panes (b) and (c) of Figs. 7 and8). Finally, the real parts of unstable eigenvalues are increased too in configuration 2 which implies that instability develops itself quicker. These sets of results tend to demonstrate that, for the studied system, a more stiff contact between the pad and the disc is more likely to lead to a stable brake behavior but that unstable ranges are harder to handle as a greater number of modes can be simultaneously unstable and hence optimization should take into account more constraints; moreover instability growth rate can become greater (larger positive real parts for eigenvalues).

The last part of this section is dedicated to the efficiency of the method compared to Monte Carlo simulations and the efficiency of each version of this method compared to the other. To this end, Table 3 summarizes CPU times for full simulations over a sample with size 1000 and 10000. For each of the proposed method version, total duration t total is split into t fit , that is the time needed for getting the decomposition and t shoot which is the amount of CPU seconds spent for evaluating the 1000 or 10000 eigenmodes using their decomposition. This table shows that the tested method is not advantageous if used on a small sample only (1000 points) for a single input law but may become so if different laws are to be tested. However as soon as a large sample (10000 points) is considered, the method is twice to four times faster than Monte Carlo simulations for a single sample, this ratio increasing with the number of samples or distributions tested. Being also accurate as demonstrated in Section 3.2, this makes it a very interesting method in terms of computational resources consumption. This table also points out a large advantage for the second version of the method based on continuation principles which is almost twice faster than the first version based on the more simple dichotomy approach.

Conclusions

This study takes place in the context of uncertain systems. Here, brake squeal noise is studied, viewed as an instability of some modes of the linear system and hence driven by the underlying eigenvalue problem. Uncertain systems can be hard to handle because of the the computation duration in the case of Monte Carlo simulations; hence a class of methods were developed to attempt to avoid the problem of repeated resolutions of a same problem with different values for uncertain parameters. A polynomial basis whose polynomials are functions of the random variables, is used to decompose the eigenmodes -eigenvalues and eigenvectors: this is the Polynomial Chaos, introduced by Ghanem and Spanos in 1991 [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] in the field of structural dynamics. Having the quantities of interest described continuously via polynomials, it is then easy to evaluate them over large samples and hence to compute the statistical estimators such as the mean or higher order moments. Depending on the problem, one can be confronted with quantities whose variations with respect to the random parameters are not continuous, that is all the more true when studying changes of stability of a dynamic system. Hence, the degree of the polynomial basis must be very large. The size of the system to solve growing with the maximum degree of the chosen basis, one can switch to another strategy and divide the stochastic space, then using a lower degree on each element of the partition: this Multi-Element generalized Polynomial Chaos (MEgPC) method was recently introduced by previous papers, mostly to study equilibria or limit cycles via time integration methods. The present work applied this approach to a structural dynamic system -pad and disc assembly -whose stability is investigated when either the friction coefficient or the contact stiffness is uncertain: the quantities of interest are the eigenmodes, solutions of the eigenvalue problem. The MEgPC method introduced in fluid dynamics is modified in two ways: the classical dichotomy approach to create the partition of the stochastic space is compared to a new version of the method based on continuation principles which happens to be more efficient. Moreover, a criterion to check the decomposition accuracy over the processed cell is proposed. This criterion, developed specifically for the stochastic eigenvalue problem is based on the Rayleigh ratio and lets us evaluate accuracy of the decomposition over a cell without any comparison to the exact values, which makes it an efficient criterion.

Both versions of the method were applied to a reduced finite element model of a brake, consisting of a pad in contact with a disc. Variations of two important yet hard to characterize interface parameters were investigated for two nominal configurations: friction coefficient and contact stiffness; results obtained via the proposed method were compared to Monte Carlo simulations. As the developed method returns a decomposition that does not depend on the random input law, different distributions were tested once the decomposition was obtained. This highlighted both accuracy and efficiency of the method, especially the version implementing continuation principles.

Post-processing of the results are to be examined in two ways. First, the choice of a random law for description of the varying parameter impacts greatly the stability profile of the brake; as characterizing this random entry is not always easy to perform, the proposed tool presents the advantage to enhance the possibility to test different input laws in a few seconds. Furthermore, the random law used can be whatever is physically thought relevant with no restriction such as being a standard normal law. Second, comparing the different contact stiffness results, one can conclude that, for the current system, the more stiff the contact is, the less the brake exhibits a squeal propensity but the more unstable modes can coexist outside of stable ranges and the greater the instability growth rate is.

Figure 1 :

 1 Figure 1: Limitations of a simple polynomial approximation of a target curve: target curve; polynomial approximation y = ψ(ξ) where ψ has degree D. (a) D = 2, (b) D = 4, (c) D = 10.

InitializationFigure 2 :Figure 3 :

 23 Figure 2: Global algorithm flowchart.Steps implying different processes for the method relying on an dichotomy approach and the one taking advantage of continuation principles are starred.

Figure 4 :

 4 Figure 4: Brake model. (a) Plan; (b) FE model.

Figure 5 :

 5 Figure 5: Modal modification due to the coupling. (a) Absolute difference between MAC matrices; (b) Eigenfrequencies of the coupled system: dots denote eigenfrequencies for eigenmodes whose MAC columns maximum absolute differences are greater than 0.5.

Figure 6 :

 6 Figure 6: Unstable modes: Config 2 (k c = 2 • 10 7 N.m -1 ), Case 2 (k c varies); Variation of complex eigenvalue λ with random variable ξ for (a) mode 13, (b) mode 80, (c) mode 82, (d) mode 85, (e) mode 86; stable range, unstable range.

  a.1) for configuration 1 and Fig. 10 (a.2) for configuration 2. For each case and each law, the following results are given for configuration 1 (k c = 1 • 10 7 N.m -1 , pane (b) of Figs. 7 to 10) and 2 (k c = 2 • 10 7 N.m -1 , pane (c) of Figs. 7 to 10

Figure 7 :

 7 Figure 7: Simulations using a uniform-like law for case 1 (µ varies): histograms. (a) Random parameter value occurrences (%); System stability and squeal information for (b.n) configuration 1 and (c.n) configuration 2; (b/c.1) Histogram of stable/unstable occurrences (%), (b/c.2) Histogram of unstable frequencies (kHz) occurrences (%), (b/c.3) Histogram of real part of unstable modes occurrences (%).

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Simulations using a normal law for case 1 (µ varies): histograms. (a) Random parameter value occurrences (%); System stability and squeal information for (b.n) configuration 1 and (c.n) configuration 2; (b/c.1) Histogram of stable/unstable occurrences (%), (b/c.2) Histogram of unstable frequencies (kHz) occurrences (%), (b/c.3) Histogram of real part of unstable modes occurrences (%).

Table 1 :

 1 Unstable ranges per mode when µ varies for (a) configuration 1 (k c = 1 • 10 7 N.m -1 ) and (b) configuration 2 (k c = 2 • 10 7 N.m -1 ).

	(a)
	Mode	µ	Im(λ) (rad.s -1 )
	6 12 69 78 81	[0.30 , µ max ] [µ min , 0.16] [9.21 • 10 3 , 9.22 • 10 3 ] [5.59 • 10 3 , 5.65 • 10 3 ] [0.20 , µ max ] [9.25 • 10 3 , 9.48 • 10 3 ] [µ min , 0.13] [2.47 • 10 5 , 2.47 • 10 5 ] [0.26 , µ max ] [3.54 • 10 5 , 3.57 • 10 5 ] [0.33 , µ max ] [3.78 • 10 5 , 3.79 • 10 5 ]
	(b)
	Mode	µ	Im(λ) (rad.s -1 )
	6	[0.29 , µ

max ] [6.10 • 10 3 , 6.28 • 10 3 ] 9 [0.31 , µ max ] [6.41 • 10 3 , 6.74 • 10 3 ] 13 [0.29 , µ max ] [9.92 • 10 3 , 1.01 • 10 4 ] 69 [0.35 , µ max ] [2.63 • 10 5 , 2.63 • 10 5 ] 73 [0.35 , µ max ] [2.92 • 10 5 , 2.92 • 10 5 ] 82 [0.20 , µ max ] [4.06 • 10 5 , 4.06 • 10 5 ] 86 [0.27 , µ max ] [4.75 • 10 5 , 4.76 • 10 5 ]

Table 2 :

 2 Accuracy of the method: Worst value over the unstable modes set for each configuration and case considering (a) method based on dichotomy and (b) method based on continuation principles.

	(a)
		Configuration 1	Configuration 2
	Indicator	Case 1	Case 2	Case 1	Case 2
	%(e Reλ < 2%)	97	76	91	75
	%(e Reλ < 5%)	98	97	97	97
	%(e Imλ < 2%)	100	87	100	97
	%(e Imλ < 5%)	100	100	100	100
	%(e φ < 5%)	100	78	100	78
	%(e sign(Reλ) )	0	0.5	0.5	0.5
	e ξHopf	0.00062 0.0097 0.0042 0.0059
	(b)
		Configuration 1	Configuration 2
	Indicator	Case 1	Case 2	Case 1	Case 2
	%(e Reλ < 2%)	97	94	91	97
	%(e Reλ < 5%)	97	97	97	99
	%(e Imλ < 2%)	100	100	100	100
	%(e Imλ < 5%)	100	100	100	100
	%(e φ < 5%)	100	100	100	100
	%(e sign(Reλ) )	0	0.5	0.5	0.5
	e ξHopf	0.00092 0.0097 0.0045 0.0012

Table 3 :

 3 CPU times comparison between Monte Carlo (MC) simulations and both versions of the proposed method: (a) method based on dichotomy and (b) method based on continuation principles; t fit : time for getting the decomposition, t shoot : time for evaluating the occurrences.

	Number	MC	(a) Dichotomy (s)	(b) Continuation (s)
	of points	(s)	t fit + t shoot = t total	t fit + t shoot = t total
	1000	345	1455 + 4 = 1459	788 + 5 = 793
	10000	3451 1455 + 44 = 1498	788 + 46 = 834

Acknowledgments

The author gratefully acknowledge the financial support of the French National Research Agency through the Young Researcher program ANR-07-JCJC-0059-01-CSD2.