Florian Hermet

Jérémie Gressier

Nicolas Binder

Shock motion inside a varying cross-section channel and consequences on the downstream flow

Keywords: Shock wave, varying cross section, A-M relation, shock decay, established flow model

 through the derivation of a one-dimensional approach connecting the shock Mach number and the area channel: the A-M relation. It is based on strong theoretical restrictions: (i) shock equations applied on a C + , (ii) omission of the post-shock influence, and (iii) initial conditions at rest. It has been the focus of many studies aimed at generalizing it. However, very little attention has been paid to the study of the shock motion outside the varying cross-sectional region since Russell (1967). The objective of the current work is to describe and to explain the shock wave behaviour in a constant area channel behind a convergent or divergent channel. It is found that the shock propagation in the downstream uniform area region is influenced by the post-shock flow unsteadinesses. Thence, the Whitham model is not suited to the shock motion study in a constant area region downstream of a convergent or divergent area region. A detailed flow description is provided and a quasi-steady model for determining the waves intensity at large times is proposed. This model gives accurate results without any assumptions on the shock strength, area variation or the shock upstream state. Finally, this study points out the limits of the use of Whitham's theory in a variable area channel with increasing section variation rate (where |dA/A(x + dx)| ≥ |dA/A(x)|).

I. INTRODUCTION

Shock wave propagation in a variable cross-section channel is a recurrent issue in the literature. Many practical occurrence of such an academic configuration exist since it happens in periodic or transient regimes of applications involving compressible gas in ducts. The flow inside a hypervelocity shock tunnel [START_REF] Igra | Simulation of the starting flow in a wedge-like nozzle[END_REF] is a typical example, but many others are found in the industry. For example in the field of aeronautical propulsion, recent discussions about the evolution of the turbofan architecture, possibly toward Humphrey's thermodynamic cycle [START_REF] Roy | Pulse detonation propulsion: challenges, current status, and future perspective[END_REF], have stimulated the interest of the scientific community concerned by turbomachinery flows for sharp transient regimes. Those would occur as a consequence of an isochoric combustion upstream of the turbines stages, which in turn would be subject to unsteady feeding.

The corresponding prototypical flow consists in the propagation of a shock-wave inside a convergent or divergent section, connecting two straight channels of different heights or diameters. The initial state of the compressible fluid can be at rest, or in a steady-state regime. The first mention of such a configuration in the literature is found in the seminal works of Chester [START_REF] Chester | The quasi-cylindrical shock tube[END_REF] and Chisnell [START_REF] Chisnell | The Normal Motion of a Shock Wave through a Non-Uniform One-Dimensional Medium[END_REF], in which the derivation of the Euler equations reveals the generic behavior expected for the shock-wave. The shock Mach number should increase inside a convergent channel (and decrease in a divergent one), which in turn produces an entropy gradient and a complex system of reflected waves. This has been experimentally confirmed * florian.hermet@isae-supaero.fr a few years later by Bird [START_REF] Bird | The effect of wall shape on the degree of reinforcement of a shock wave moving into a converging channel[END_REF]. On the modelling ground, the early proposition of Whitham [START_REF] Whitham | On the propagation of shock waves through regions of non-uniform area or flow[END_REF] still holds a lot of attention in the scientific community. The theoretical problem is derived through a one-dimensional approach of the shock wave propagation in a gas at rest inside a variable cross-section channel. The entry point of this derivation is the solving of the shock equations along a forward propagating characteristic (C + ). It formalizes the so-called A-M relation, which relates the shock Mach number, M s , to the local cross-section area of the channel, A. The A-M relation thus predicts the evolution of the shock intensity during its propagation in a convergent or a divergent channel. The first derivation of this relation was restricted to smooth variations of the cross-section. An extension for more abrupt geometries is possible, if the effects of the post-shock flow on shock motion are neglected (Whitham [8]).

The A-M relation thus provides a comprehensive modelling of the shock-wave propagation in non-uniform duct flows, widely employed in the literature. Its validity, addressed through a comparison with experimental or numerical data, is questioned in many references: the validity of the Whitham model has been discussed by Igra et al. [START_REF] Igra | Shock wave interaction with area changes in ducts[END_REF] as well as its ability to cope with additional effects such as the influence of the wall shape of the converging [START_REF] Bird | The effect of wall shape on the degree of reinforcement of a shock wave moving into a converging channel[END_REF] or diverging [START_REF] Nettleton | Shock attenuation in a gradual area expansion[END_REF][START_REF] Gottlieb | Interaction of rarefaction waves with area reductions in ducts[END_REF] channels, or even the maximization of the shock strength [START_REF] Dowse | Area change effects on shock wave propagation[END_REF].

A general consensus of the scientific community is agreed on the fair accuracy of the A-M relation for moderate shock intensity and moderate area ratio, or at least smooth evolution of the convergent channel. The limitations of the model are likely to be a consequence of the strong theoretical restrictions imposed in its the derivation: (i) shock equations applied on a C + , (ii) omission of the post-shock influence, and (iii) initial conditions at rest. In addition, it appears that a singularity exists at sonic conditions downstream of the shock wave, which has been solved by Friedman [START_REF] Friedman | An improved perturbation theory for shock waves propagating through non-uniform regions[END_REF]. Later, Yousaf [START_REF] Yousaf | The effect of overtaking disturbances on the motion of converging shock waves[END_REF] and Milton [START_REF] Milton | Mach reflection using ray-shock theory[END_REF] have derived an alternative formulation of the A-M relation taking into account the overtaking disturbance on the shock motion, using a set of assumptions. A more general theory has emerged a few decades ago and can be found in Best [START_REF] Best | A generalisation of the theory of geometrical shock dynamics[END_REF]. The author derived a system of differential equations that does not require the omission of the post-shock flow effect on the shock motion. Best model predicts the shock motion at an higher order than the A-M relation. The possibility to prescribe a different initial condition than a quiescent uniform gas is examined in Chisnell [START_REF] Chisnell | A note on whitham's rule[END_REF] and Whitham [START_REF] Whitham | A note on shock dynamics relative to a moving frame[END_REF], in which the A-M relation is adapted to the propagation of a shock wave in a moving fluid. Similarly Catherasoo and Sturtevant [START_REF] Catherasoo | Shock dynamics in non-uniform media[END_REF] revised the A-M relation for the case of a nonuniform quiescent gas.

All these studies are focused on the unsteady flow observed inside the convergent or divergent channel. It is indeed in that section that the main interaction between the shock wave and the geometry takes place. However, very little attention has been paid to the shock motion in the straight channel downstream of this critical element. A return to a normal behaviour, i.e. a shock wave with a constant intensity, would stand as a reasonable expectation since it propagates in a simple duct. However this is not the case, and this is seldom mentioned in the literature. The work of Russell [START_REF] Russell | Shock-wave strengthening by area convergence[END_REF] seems to be the only one to propose an explicit description of the shock propagation beyond the convergent or divergent section, and it is moreover treated as a side effect. In this experimental study, the author tries to quantify the strengthening of the shock intensity within the convergent channel. For this, pairs of thin-film gauges are used to determine the average speed of the shock wave entering the convergent. Three measurement positions downstream of the convergent channel are also recorded. The shock wave strengthening was actually observed to occur inside the convergent channel, and a progressive decay of the shock intensity was observed in the downstream duct. This decay is briefly analyzed in the work of Russel, and the downstream deceleration was imputed to the viscous dissipation, as well as higher order interactions.

The objective of the present work is to describe and explain the shock wave behaviour in an uniform area channel behind a convergent or divergent channel, in an inviscid frame so as to discard viscous effects. Some propositions are made in order to assess the validity of the A-M relation since it is not suitable to predict the shock wave behaviour in the downstream straight channel, i.e. at larger time. An accurate description of the flow and the waves interactions during shock wave motion is proposed. A model predicting the flow configuration at larger times is demonstrated to determine the wave intensities inside the transient flow. The proposed model does not impose any restriction in terms of initial condition of the flow (rest or moving), intensity of the initial shock, or aspect ratio of the convergent/divergent channel.

The paper is organized as follows. In the first section, we briefly review the model of Whitham before discussing its validity in the next section by using numerical simulations. The last part of the paper is dedicated to the description of the long-time flow physics, together with the validation of our model. The paper ends with the conclusions to this study.

II. REVIEW OF WHITHAM MODEL A. Physics description

We first briefly recall the main flow features when a shock wave propagates through a convergent and/or a divergent channel. We consider a varying cross-section channel which streamwise evolution is given by:

A = A in if x < 0 A = A(x) if x > 0 ( 1
)
The flow is initially at rest. A shock wave is imposed initially (t = 0) at the inlet of the channel, through the prescription of the inlet boundary conditions. In the first part of the channel (negative values of x), this shock wave propagates in a uniform quiescent gas in a constant section duct. Its propagation speed is constant and measured by the shock Mach number M s . The flow is thus divided into two uniform regions separated by the shock wave. The undisturbed state ahead of the shock is denoted by (0) while the uniform state prescribed at t = 0 behind the shock is denoted by [START_REF] Whitham | On the propagation of shock waves through regions of non-uniform area or flow[END_REF].

When the shock reaches the convergent or divergent part of the channel, its intensity is modified in response to the area variation. Therefore, disturbances induced by this modification of the shock intensity travel upstream, following negative characteristic lines (C -) and the particles paths (P ), as illustrated in figure 1 for the case of a shock wave propagating through a convergent duct.

The shock wave accelerates, and gains in intensity. Consequently, a reflected compression wave and a positive entropy gradient are generated. On the (x,t) diagram of figure 1, the flow behind the shock (1) remains subsonic. For a supersonic downstream flow, the slope of the negative characteristic C -would be positive. Conversely for a divergent duct, a reflected expansion wave is propagated by the C -and a negative entropy gradient propagates on the particles path.

As stated in the introduction, the A-M relation gives an accurate prediction of the physics of such a flow configuration. In the next subsection, we recall the theoretical derivation of this model proposed by Whitham [START_REF] Whitham | On the propagation of shock waves through regions of non-uniform area or flow[END_REF].

FIG. 1. Spatio-temporal (x,t) diagram for a shock propagating in a convergent duct (dA < 0). The flow is subsonic downstream of the shock.

B. Demonstration of the A-M relation

The demonstration proposed by Whitham is based upon the quasi one-dimensional and inviscid governing equations, written in the characteristic form downstream of the shock:

dp ρa + du + ua A dA(x) dx = 0 for C + : dx dt = u + a (2a) dp ρa -du + ua A dA(x) dx = 0 for C -: dx dt = u -a (2b) ds = 0 for P : dx dt = u (2c)
where :

d = ∂ ∂t + dx dt ∂ ∂x (3)
stands for the material derivative; p = p(M s ), ρ = ρ(M s ), u = u(M s ) and a = a(M s ) are the conditions behind the shock, given by the Rankine-Hugoniot jump relationships. Whitham suggested that the shock path can be approached by that of converging characteristics on the shock. As illustrated in the sketch of figure 1, the shock path is approximated by a forward characteristic C + . The A-M relation is formalized by the substitution of the differentiated Rankine-Hugoniot jump relations on the differential invariant of the C + , equation (2a):

1 A dA dM s = -g(M s ) (4)
with :

             g(M s ) = M s M 2 s -1 1 + 2 γ + 1 1 -µ 2 µ 1 + 2µ + 1 M 2 s µ 2 = (γ -1)M 2 s + 2 2γM 2 s -(γ -1) (5)
The validity of this strong assumption made by Whitham to approximate the shock path by a C + can be evaluated analytically. Applying the differential invariant of the C + to the shock moving at the speed a 0 M s , Whitham claims that equation ( 6) is a good approximation of the shock behaviour:

dp ρa + du + ua A dA(x) dx = 0 for C + : dx dt = a 0 M s (6)
The difference between equations (2a) and ( 6) provides a quantification of the error caused by the two assumptions made by Whitham:

1 u + a (u + a) -a 0 M s a 0 M s 1 ∂p ∂t + ρa ∂u ∂t 2 (7)
The first term gives a measure of the coincidence between the shock and the C + . For a weak shock, it is close to zero as a 0 M s ≈ u+a, however this term tends to 0.274 as M → ∞ for γ = 1.4. The second term gives a quantification of unsteadiness in the post-shock flow, which modifies the shock strength. It is purely neglected in Whitham's demonstration, and very few comments about it are found in the literature. The discussions are essentially focused on the coincidence between the shock and the C + , an assumption that gives to Whitham approach a better validity for the case of a weak shock rather than for a strong shock. However, the unsteadiness of the postshock flow can induce some errors, even for a weak shock. In the next section, we will show that the post-shock unsteadiness must be taken into account to correctly predict the shock intensity downstream of a convergent or divergent duct.

III. VALIDITY OF THE WHITHAM MODEL

According to the A-M relation, the shock motion only depends on both the initial shock strength and the area streamwise variation. A numerical integration of the onedimensional Euler equations is performed to assess the accuracy of Whitham model. The simulations are carried out thanks to in-house IC 3 solver, forked from CharLES X solver [START_REF] Bres | Unstructured large-eddy simulations of supersonic jets[END_REF] in which the area variations are modelled by volumic source terms, equation (8). IC3 is based on the resolution of the compressible formulation of Euler equations in their conservative form, spatially filtered, on an unstructured mesh using a finite volume method. An explicit third-order Runge-Kutta scheme is used for time advancement while an essentially non-oscillatory (ENO) second-order shock-capturing scheme is applied to compute flux.

               ∂ρ ∂t + ∂ρu ∂x = - ρu A dA dx ∂ρu ∂t + ∂ ρu 2 + p ∂x = - ρu 2 A dA dx ∂ρE ∂t + ∂ (u (ρE + p)) ∂x = - u (ρE + p) A dA dx (8)
The shock motion inside a duct with cross-sectional variation is thus considered for different initial strengths and several area variations. The study is performed for a linear convergent channel and restricted to cases for which the flow is subsonic behind the shock. The flow ahead of the shock is at rest, in accordance with Whitham theory.

A. Shock strength prediction inside a variable section channel

Figure 2 shows a comparative study between the prediction given by the A-M relation of Whitham and the results of our numerical simulations. Both weak shock ( P1 P0 = 1.2) and a strong shock ( P1 P0 = 3.0) configurations are compared for a small ( Ain Aout = 1.2) and an high ( Ain Aout = 3.0) converging ratio. A good agreement between the theoretical predictions and the numerical simulation can be observed for all cases. As expected, the error increases with the shock strength. The error between simulation and the A-M relation is computed as: Error = P simulation -PA-M PA-M-P0 where (0) stands for the undisturbed state ahead of the shock. The relative error increases with the area ratio. This can be explained by significant flow unsteadinesses. Those post-shock flow unsteadinesses are a consequence of the shock strengthening in the converging channel, and can be interpreted as the duct response to the brutal evolution of the flow. The unsteadinesses intensity is stronger when the shock modification is considerable, which happens for large converging ratios. The present results show good agreement with the literature and do not question the accuracy of the A-M relation for these specific flow configurations. However, very few studies have focused on what happens in the straight channel just downstream of the convergent/divergent section, in which the A-M relations predicts a constant evolution of the shock. This point is discussed in the following subsection.

B. Shock wave propagating downstream of the convergent channel

The shock motion is investigated for a single initial pressure ratio ( P1 P0 = 1.2) and single convergent channel linear ratio ( Ain Aout = 3), but now the region of interest is extended to the constant area channel following the convergent section.

Figure 3-a presents the evolution of the shock pressure ratio along its propagation inside the channel, predicted by the numerical simulation and the A-M relation. It is clear that the shock strength does not remain constant in the straight channel downstream of the convergent channel, as predicted by the A-M relation. The shock wave intensity actually decreases, since it is weakened by a forward propagating expansion wave. This expansion wave comes from the interaction between the backward propagating compression wave (the reflected wave), the entropy gradient created by the shock strengthening inside the convergent duct and convected by the flow, and the cross-section evolution. Since the backward traveling reflected wave is a compression wave that evolves in a divergent section (because of the backward direction), an expansion wave is created, this time traveling in the flow direction. All theses effects compose the complex and unsteady post-shock flow that activates the term 2 in equation ( 7).

However, it may be surprising at first sight that the post-shock flow effects have more influence on the shock motion in the downstream straight channel, rather than in the convergent/divergent section. Actually it is not as explained by Best [START_REF] Best | A generalisation of the theory of geometrical shock dynamics[END_REF].

Best demonstrated a valid recurrence relation along the shock path predicting the shock wave intensity. This relation is a handling of the Euler equations, which does not require any of the aforementioned assumptions used by Whitham; for more details see Best [START_REF] Best | A generalisation of the theory of geometrical shock dynamics[END_REF]. The final equations of Best model are:

dM s = -a 0 M s d Ms p + ρad Ms u ρa 2 u u + a 1 A dA dx + 1 u + a - 1 a 0 M s Q 1 (9a) ∀k ≥ 1, dQ k = -a 0 M s ∂ k ∂t k ρa 2 u u + a 1 A dA dx + k i=1 k i ∂ i ∂t i 1 u + a Q k-i+1 (9b) + ∂ k-1 ∂t k-1 ∂ρa ∂t ∂u ∂x - ∂ρa ∂x ∂u ∂t - a 0 M s u + a -1 Q k+1
                 d = ∂ ∂t + a 0 M s ∂ ∂x Q k+1 = ∂ k ∂t k ∂p ∂t + ρa ∂u ∂t , ∀k ≥ 0 d Ms p (resp. u) represents a derivative of p(M s ) (resp. u(M s )) to M s , (10)
and show that the truncation of Best's recurrence system to the first equation (i.e Q 1 = 0 and consequently dQ k = 0) allows to get back exactly the A-M relation proposed by Whitham. Thus, a reliability criterion of the A-M relation can be given, based on the equation (9a):

a 0 M s ρa 2 u u + a 1 A dA dx ACE a 0 M s u + a -1 ∂p ∂t + ρa ∂u ∂t PSE (11)
The left hand side of equation ( 11) depicts the area change effects (ACE) on the shock strength. The right hand side has already been described in the previous section as the error carried out by the A-M relation, and corresponds to the post shock effects (PSE). From equation [START_REF] Gottlieb | Interaction of rarefaction waves with area reductions in ducts[END_REF], it is clear that the Whitham model is accurate when the area change effects are predominant over the shock propagation, compared to the post-shock flow effects.

The post-processing of the numerical simulations gives access to the exact balance between the two physical mechanisms: the area change effect (ACE) and the unsteady post-shock effect (PSE). Figure 3-b shows the evolution of each of the two terms, immediately behind the shock; following its propagation. Obviously, ACE = 0 for each straight part of the channel, which actually makes those regions very sensitive to post-shock effects. In the first straight part of the channel, the shock propagation is uniform, which makes PSE = 0 and, as a consequence, the Whitham model is exact. In the converging section PSE is activated, but so does ACE, which is one or even two orders of magnitude higher than PSE: the error of the model exists but it is not perceptible. Consequently, Whitham theory fits very well the simulations in that part of the channel. Then, in the downstream straight section, PSE is still activated because of the expansion wave traveling upstream, but ACE = 0: the error gets significant and the Whitham model is inaccurate. Finally, effects due to the post-shock flow disturbances tend towards zero as the shock strength converges to its asymptotic value in the uniform channel.

The truncation at the 2 nd equation of Best system is also examined in figure 3. In the converging section, no real improvement is to be expected. In the downstream straight channel, a decay of the shock intensity is effectively predicted, since the post-shock effects are partially computed. But the truncation makes the prediction at long time (i.e. when the interaction between the shock and the expansion wave is over) quite inaccurate, compared with the direct simulation. An higher order truncation or a closure model would allow getting a more pertinent strength at long time by better modelling the post-shock unsteadinesses (equation (9a)). However, the recurrence system demonstrated by Best (equation ( 9)) is quite complex to solve at a higher order: to the authors knowledge, no implementation of Best model with k ≥ 2 or with a closure model has ever been carried out.

For completeness, it must be stated that some models of the literature, such as found in Yousaf [START_REF] Yousaf | The effect of overtaking disturbances on the motion of converging shock waves[END_REF] and Milton [START_REF] Milton | Mach reflection using ray-shock theory[END_REF], also include the post-shock flow effects on the shock motion, but this modeling is only active in the converging/diverging section (when dA = 0 the shock strength no longer varies in their modelling). This is unfortunate since no correction of the Whitham model is needed in that region, as opposed to the downstream constant section channel of this region for which the post-shock flow effects should be taken into account in order to capture the shock decay.

C. Non-linear area variation in the converging/diverging section

In the previous simulations, A(x) is linear. This provokes a brutal cancellation of the ACE term at the inlet of the straight channel, which might not be representative of more realistic cases, where the geometric interface might be more progressive. Some cases for which dA varies in the converging section are examined in that purpose. The initial strength of the shock and the global area variation are the same as previously P1 P0 = 1.2 and Ain Aout = 3.0 . Figure 4 shows the influence of the sign of d 2 A on the intensity of the ACE term. Actually, the result observed does not drastically change the former conclusions: the error gets significant when the condition given by equation [START_REF] Gottlieb | Interaction of rarefaction waves with area reductions in ducts[END_REF] is no more respected. Anyway, it moderates the assessment that the Whitham model is valid inside the converging section since, technically, this condition is violated while the shock wave is still inside the convergent channel. A strong correlation is found between the variation of the ACE and that of d

IV. PHYSICS DESCRIPTION

The physical framework of the model is recalled, illustrated by the specific case of a shock wave propagating inside a linear convergent channel followed by a straight channel; the flow downstream the shock is still considered as subsonic. But this focus on a given configuration does not alter the generic nature of the model, as will be shown in the following sections.

The physical description proposed in section III B is illustrated in the (x,t) diagram of figure 5. The key features are briefly summarized:

• the shock strengthens in a convergent (resp. weakens in a divergent) channel due to the cross-section change;

• a compression wave is reflected (resp. an expansion wave for dA < 0) in the post-shock flow, travelling at u -a;

• an entropy gradient is created by the varying shock intensity, travelling at post-shock flow speed u;

• in the post-shock flow, the interaction between the reflected wave and the entropy gradient and their respective interaction with the area change, produce the disturbances that activate the PSE term in equation [START_REF] Gottlieb | Interaction of rarefaction waves with area reductions in ducts[END_REF].

FIG. 5. The (x,t) diagram of the shock wave propagation through a linear convergent channel followed by a straight channel

Inside the convergent section, the ACE term is more than an order of magnitude higher than the PSE term, at least as long as |dA(x + dx)| ≥ |dA(x)|, which is the case for a linear area variation. However, these disturbances weaken the shock wave in the downstream uniform crosssection region, which is precisely what needs to be modeled.

Figure 5 details the post-shock flow unsteadiness, which generates the expansion wave propagation along the C + lines that finally hits the shock. The main reason for the appearance of this expansion wave is the interaction between the reflected wave and the area change: for a convergent duct, the reflected compression wave interacts with a divergent channel. Thus, the compression wave weakens and induces a streamwise expansion wave. The exact influence of the entropy gradient is less clear, but this entropy gradient also persists in the straight channel, since the shock wave intensity is decaying. Regarding this decay, it is worth noting that it ends after a characteristic length that is not only correlated to the area ratio Ain

Aout , but also to the axial length of the convergent/divergent section and the shock wave intensity. Moreover, no reflected wave is considered by the shock decay, as mentioned by Courant and Friedrichs [START_REF] Courant | Supersonic Flow and Shock Waves[END_REF].

Finally, the disturbance carried by the C + (here the expansion wave) interacts in turn with the section variation. This produces a new reflected wave that will follow the initial one and will, in turn, interacts with the area change until the total attenuation of disturbances carried by C + and C -in the convergent/divergent section is reached. This equilibrium state requires many iterations for very intense cross-sectional variations Ain Aout > 10 . For less severe area ratios, the observations show that the first reflection of the reflected wave (here an expansion), carried by C + , almost directly sets the gas inside the convergent channel in a state compatible with the quasisteady approach (Nozzle theory, initiated by De Laval [START_REF] Laval | Steam turbine[END_REF]). Figure 6 illustrates this observation. The time t 1 , also quoted in figure 5, corresponds to the position of the shock wave exactly at the outlet of the converging section. The state of the gas inside the convergent duct is not compatible with the quasi-steady approach, and the same is true for a moment before t 2 . At t 2 , which corresponds to an instant following the departure of the last C + of the streamwise expansion wave from the convergent channel, the flow is adapted to the quasi-steady solution.

The complete attenuation of the wave interactions inside the convergent channel, together with a stabilization of the shock wave intensity, marks the beginning of the long time state (t 3 , on figure 5 &6). At long time, states (2), ( 3) and (4) (figure 5) persist in the flow (figure 6), until a possible wave interaction with the boundaries of the straight channel. The next section is dedicated to the determination of these states.

V. ESTABLISHED FLOW MODEL

A. Model description

The model is based on the aforementioned observation: the flow in the convergent/divergent section can be approximated by the quasi-steady approach as soon as the streamwise expansion leaves this section. The validity should thus be restricted to area-ratios lower than 10, but the following section will show that it provides fair prediction even above that limit.

The states (1) and ( 2) are separated by an isentropic compression wave (figure 5) in a straight section, so the Riemann invariant along a C + characteristic for a perfect gas together with the isentropic relation gives:

         u 1 + 2a 1 γ -1 = u 2 + 2a 2 γ -1 P 1 a 2γ γ-1 1 = P 2 a 2γ γ-1 2 (12)
The flow inside the convergent channel is here considered as quasi-steady, as formerly stated. Therefore, the classical steady relations of mass flow and entropy conservation can be applied between the states (2) and ( 3), where P i represents the stagnation pressure.

                       A in A out = M 3 1 + γ-1 2 M 2 2 γ+1 2(γ-1) M 2 1 + γ-1 2 M 2 3 γ+1 2(γ-1) P 3 = P i2 1 + γ -1 2 M 2 3 -1 γ-1 P 2 a 2γ γ-1 2 = P 3 a 2γ γ-1 3 (13)
The states (3) and ( 4) are only separated by an entropy gradient. On each side of this gradient, conditions are the same as those of a contact discontinuity, hence:

P 4 = P 3 and u 4 = u 3 (14)
Thanks to the Rankine-Hugoniot equations, which can be applied between states (0) and (4), state ( 4) is determined by shock intensity or velocity. The right number of relations has been obtained to close the problem and determine states (2), (3) and (4). Thus, all long-time states can be determined through an iterative process:

• initialization of the reflected wave intensity between (1) and (2), state (2) is then fully determined;

• nozzle steady relations provide state (3);

• contact discontinuity relations enforce velocity and pressure of state (4); FIG. 6. Density is plotted for the specified times on Fig. 5. P 1 P 0 = 1.2, A in A out = 3.0. Red, green and blue circles are respectively the nozzle steady relation computed from the final steady state, instantaneous states and the state just behind the shock wave, in the outlet cross-section of the convergent channel • shock intensity is then determined to match the pressure condition of state (4) through Rankine-Hugoniot relations, which also prescribe the entropy gradient and the velocity;

• the post-shock velocity of state (4) must match the already prescribed velocity; otherwise, the initial guess for the reflected wave strength is wrong and must be modified accordingly until convergence.

No assumptions about the shock initial strength, the cross-sectional variation intensity or the upstream state to the shock are required. Furthermore, the maximum value of entropy in the flow, between the positive and negative entropy gradients, can be estimated using the A-M relation or Best model for a quiescent gas ahead of the shock. For a moving flow, it can be estimated with Chisnell's model [START_REF] Chisnell | A note on whitham's rule[END_REF].

B. Model validation

In order to validate the model, comparisons with numerical simulations are performed inside a linear convergent channel located between two straight sections, for several initial shock intensities. Comparisons are achieved by examining the density, since this quantity varies across states (2), ( 3) and (4) (figure 5) unlike the pressure or velocity, for which a continuity is expected across the entropy gradient.

For a fixed area-ratio and several initial shock strengths (graphs on the left in figure 7), the model gives very good results. The long-time transmitted shock wave is very well predicted since the error is less than 1%. There is no clear trend regarding the evolution of the error between the numerical simulation and the The model shows also excellent results for a range of converging ratio at constant initial shock intensity (graphs on the right in figure 7). However, the error tends to increase for large ratio, as expected above 10. However, even for large cross-sectional variations, the re-sults are acceptable (Errors < 2%).

The validation presented here is restricted to a linear variation inside a convergent duct, for compactness. However, similar results have been obtained for nonlinear section variations, and for divergent sections. Indeed, no hypothesis is required on the section variation shape in the equation [START_REF] Friedman | An improved perturbation theory for shock waves propagating through non-uniform regions[END_REF]. And even if the physics actually change for a divergent (the reflected compression wave is replaced by a reflected expansion wave), it has no consequences regarding the model since the equations through an isentropic compression or an isentropic expansion wave remain the same. Therefore, the model still holds.

VI. CONCLUSION

This paper highlights the physics of the flow created by the shock wave propagation through a variable area channel. Current models are assessed by resolved numerical simulations and their limitations are demonstrated, more specifically in the straight section downstream of a converging (or diverging) region. The paper demonstrates that the shock wave strength can be drastically modified by the post-shock flow unsteadinesses, which contrasts with Whitham's assumptions in the derivation of the well-known A-M relation. The mechanism of interaction between the reflected waves, which provoke this postshock unsteadiness, and the shock wave itself is described in the paper, and illustrated by the use of (x,t) diagrams. Its influence on the shock wave intensity is hardly perceptible inside the converging/diverging section, which makes Whitham model quite valid in that region. But it strongly increases in the downstream straight channel, and prevents any accurate prediction at long time with existing models of the literature, including Best model with two equations. A higher truncation of Best model would probably enhance its accuracy, but the cost of resolution becomes comparable to that of a direct simulation.

A simple quasi-steady model is also proposed in this paper. This model accurately predicts the long time waves intensity, without any hypothesis on initial shock wave strength or regarding the area variation intensity, as well as on the shock upstream flow state (rest or moving). The validity of the model is assessed on the specific case of a linearly converging section in this paper, but it has been observed with equal accuracy for non-linear evolution of the area and/or diverging sections.

FIG. 2 .

 2 FIG. 2. Comparison between the numerical simulation and the A-M relation for various area ratio and initial shock strength.

  2 A. However, a generic criterion based on the geometry variation seems out of the reach of analytical derivation. But it could be guessed that the Whitham model gives relatively accurate results when |dA(x + dx)| ≥ |dA(x)| and cannot provide good results when |dA(x + dx)| < |dA(x)|. Unfortunately, the condition |dA(x + dx)| ≥ |dA(x)| never holds at the end of the convergent/divergent duct, where the channel becomes straight, no matter how smoothly. A specific model is thus required to predict the decay of the shock wave, caused by PSE term. Such a long-time model is proposed in the following section.

FIG. 4 .

 4 FIG. 4. Investigation of the d 2 A variation through the non-uniform region. (a): shock strength during its propagation. (b): terms of equation (11), each term is normalized by the maximum of ACE. (c): the area evolution along the geometry.

FIG. 7 .

 7 FIG. 7. Comparison between numerical simulation and long time model for different initial shock strength and for various linear area change.

.

ACKNOWLEDGEMENTS

The authors would like to thank the DGA for their financial support and for permitting the publication of the research. This work was performed using HPC resources from GENCI-IDRIS (Grant 2020-100431), GENCI-CINES (Grant 2020-A07178) and CALMIP (Grant 2020-p1425).


Cherif Snoun
email: a.cherif.snoun@insa-cvl.fr

Baptiste Bergeot
email: b.baptiste.bergeot@insa-cvl.fr

Sébastien Berger
email: c.sebastien.berger@insa-cvl.fr

Sébastien Berger Generalized Polynomial

Generalized Polynomial Chaos for robust modelling of Nonlinear Energy Sinks used to mitigate dynamic instabilities in braking systems

Keywords: Friction-induced vibration, Nonlinear Energy Sinks, Uncertainty, Robust approach, generalized Polynomial Chaos

Introduction

In the working process of the braking system, squealing noise may be produced if the system comes into an unstable state. It may affect the user comfort and increases the customer complaints. Developing strategies to better understanding and control instabilities is therefore crucial.

Several researchers have been developed to understand the mechanisms generating noise and to predict it [START_REF] Kinkaid | Automotive disc brake squeal[END_REF][START_REF] Tison | Improvement in the predictivity of squeal simulations: Uncertainty and robustness[END_REF][START_REF] Wang | Noise performance improvements and tribological consequences of a pad-on-disc system through groove-textured disc surface[END_REF]. Friction-induced vibrations or dynamic instabilities are the scientific terms almost used. These phenomena can appear by the generation of Limit Cycle Oscillation (LCO) induced by dry friction. They are explained, in most studies of self-excited systems, by the coupling of the tangential and normal modes [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF][START_REF] Hervé | Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects[END_REF]. As a way to model these dynamic instabilities related to friction, the wellknown two degrees of freedom Hultèn's model [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF][START_REF] Hultén | Brake squeal -a self-exciting mechanism with constant friction[END_REF] has been widely used and it is also considred as the primary system in this paper.

In order to mitigate LCO, we propose to use the concept of Targeted Energy Transfer (TET). In the last decades, TET has become an important passive control technique for reducing or eliminating unwanted vibrations [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF]. It consists to attach a nonlinear device also named Nonlinear Energy Sink (NES) to the main system. The NES represents a nonlinear spring mass damper with a nonlinear stiffness, typically cubic as considered in this work. It can adapt itself to the main system that it is attached to without being tuned to a specific frequency.

NES has been widely studied theoretically and experimentally. It has been applied in several applications as for vibration mitigation [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF], noise reduction [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF] or seismic mitigation [START_REF] Nucera | Targeted energy transfers in vibro-impact oscillators for seismic mitigation[END_REF]. In recent years, particular attention has been paid for the NES design in deterministic approach [START_REF] Tuan | Design criteria for optimally tuned nonlinear energy sinks-part 1: transient regime[END_REF] and also in probabilistic approach incorporating uncertainty in the parameters of the model [START_REF] Boroson | Optimization under uncertainty of parallel nonlinear energy sinks[END_REF]. In fact, uncertainty might have a large effect on the performance of NES. In the previous work by the authors [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF], two ungrounded NESs are coupled to the well-known Hultèn's model in order to mitigate or eliminate squeaking noise in the braking system. Related to the dispersion of some uncertain parameters (here the friction coefficient), the authors classified the steady-state response regimes in two main regimes: the first is mitigated regimes which are composed by complete suppression regimes, periodic response regimes and strongly modulated responses regimes and the second is unmitigated regimes. The authors showed that the prediction of the boundary value between strongly modulated responses and no mitigation responses is not performed. For that, it becomes necessary to take into account these uncertainties in the study of the dynamical behavior of the system.

In the literature, several probabilistic approaches are used for propagating uncertainties in a deterministic model (DM) of a mechanical system and then to evaluate the quantity of Interest (QoI) which is a response from the DM. The generalized Polynomial Chaos (gPC) method is a well-known metamodel using for quantification of uncertainties. It is a less costly alternative to the Monte Carlo (MC) approach which is often used as reference [START_REF] George S Fishman | Monte carlo: concepts, algorithms, and applications[END_REF]. The gPC method expands the QoI of the DM in series of polynomials of uncertain variables. Many papers using the gPC method have been published as for the stability analysis of a break system [START_REF] Nechak | Non-intrusive generalized polynomial chaos for the robust stability analysis of uncertain nonlinear dynamic friction systems[END_REF][START_REF] Nechak | Robustness of stochastic expansions for the stability of uncertain nonlinear dynamical systems -application to brake squeal[END_REF] and more recently for the stability analysis of a clutch system [START_REF] Thinh Kieu | Stability analysis of a clutch system with uncertain parameters using sparse polynomial chaos expansions[END_REF].

In this paper, the well-known Hultèn's model coupled to two ungrounded NESs is used. As explained above, a friction system coupled to NES may be in a regime with harmful LCO or not related to the dispersion of the uncertain parameters. Moreover, the transition from the mitigated regimes to the unmitigated regimes implies a jump in the steady-state amplitude profiles. Thus, the efficiency of the NES to reduce or eliminate squealing noise can switch from high to low for a slight change of the uncertain parameters. The aim of this paper is therefore to develop a method for robust modelling NES by locating this discontinuity in the dynamical behavior of the system and then to predict, with a low computational cost, the Propensity of the PHSSR. For that, two techniques using generalized Polynomial Chaos are investigated.

The article is constructed as follows: In Section 2, the gPC theory is presented. In Section 3, the two degrees-of-freedom selfexcited model coupled to NESs is introduced. In Section 4, the problem statement is formulated. In Section 5, the two jump detection methods using the gPC are developed. In Section 6, the results of the method presented in Section 5 is shown. Finally, conclusions are given in Section 7.

generalized Polynomial Chaos theory 2.1 generalized Polynomial Chaos

The generalized Polynomial Chaos (gPC) theory has been introduced by Wiener [START_REF] Wiener | The homogeneous chaos[END_REF] and launched by Ghanem and Spanos [START_REF] Roger | Stochastic finite elements: A spectral approach[END_REF]. The representation in the base of gPC called also Wiener-Hermite development, allows a description of the random second order function. Consider (Ω,β,Pr) a probability space where Ω is the sample space of the random events ω, β is the σ-algebra of the subsets of Ω and Pr is the probability measure.

Let µ(µ 1 , ..., µ r ) is a vector of r real parameters supposed to be uniformly distributed within a given interval [a, b] r and ξ(ξ 1 , ..., ξ r ) is a vector of r independent random variables within[-1, 1] r and obtained as

µ j (ξ j ) = a j + b j 2 + b j -a j 2 ξ j , (j = 1, . . . , r). (1)
The generalized Cameron-Martin theorem [START_REF] Cameron | The orthogonal development of non-linear functionals in series of fourier-hermite functionals[END_REF] determined that it is possible to express a random process X(ξ) called also the Quantity of Interest ( QoI ) with a truncated orthogonal polynomial function series such as

X(ξ) ≈ Np j=0 xj φ j (ξ), (2)
where xj are the gPC coefficients of the stochastic process X(ξ), φ j (ξ) are orthogonal polynomial functions and N p = (p+r)! p!r! -1 where p is the order of the gPC and r is the number of the uncertain parameters [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF].

The polynomials φ j is an orthogonal basis in the Hilbert space L 2 (Ω, β, P r) with the following orthogonality relation:

φ i , φ j = φ i , φ i δ ij (3)
where δ ij denotes the Kronecker delta

δ ij =    0 if i = j 1 if i = j
and φ i , φ i = φ i φ i W (ξ)dξ where W (ξ) is the probability density function of ξ.

In this paper, the gPC coefficients are determined using a non-intrusive methods from a quite number of values of the QoI, obtained from numerical simulations of the determinstic Model (DM) presented in the Section 3. The non-intrusive method used is the regression approach which requires a minimum of Q = kN p (where k is a small integer usually equal to 2 or 3) simulations to built the coefficients. The Latin Hypercube Samples (LHS) method [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] are used to perform these Q simulations.

Using the regression approach, the coefficients are determined through the minimization of the criterion [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF] given by

ε 2 reg = Q q=1   X ξ (q) - Np j=0 xj φ j ξ (q)   2 , (4)
where

ξ (q) = ξ (q) 1 , ..., ξ (q) r
, with q = 1, ..., Q, is the Numerical Experimental Design (NED) chosen from the LHS method. X ξ (q) is the vector of the corresponding evaluations of the DM.

Finally, the gPC coefficients xj are defined by

x = φ T (ξ (q) )φ(ξ (q) ) -1 φ T (ξ (q) )X(ξ (q) ), (5)
where x = (x j , . . . , xQ ) and φ(ξ (q) ) is the matrix defined by

φ(ξ (q) ) =     φ 0 (ξ (1) ) . . . φ Np-1 (ξ (1) ) . . . . . . . . . φ 0 (ξ (Q) ) . . . φ Np-1 (ξ (Q) )     . (6)

Criterion for choosing the order of truncation

In this section, we present the criteria for choosing the optimal order of the truncated polynomial chaos expansion (see Eq. ( 2)). To obtain a good accuracy of the gPC, a high polynomial order is required and consequently a high computational cost. It is therefore necessary to determine the optimal order of truncation in terms of accuracy and computational cost.

Let xj,p and xj,p-1 the vectors of the gPC coefficients evaluating respectively by two successive order developments p and p -1. The corresponding evaluations of the QoI by the gPC are given by:

             X p (ξ) ≈ Np j=0 xj φ j (ξ) X p-1 (ξ) ≈ N p-1 j=0 xj φ j (ξ) (7)
According to the statistical characteristics of the gPC, we have:

                   xp = xp-1 = x0 σ 2 p = 1 2 r Np j=1 x2 j φ 2 j σ 2 p-1 = 1 2 r N p-1 j=1 x2 j φ 2 j ( 8
)
where xp is the mean and σ 2 p is the variance of the QoI by the gPC in the order p. The error of variance e σ 2 ,p between two successive orders gPC is given by:

e σ 2 ,p = σ 2 p -σ 2 p-1 = 1 2 r Np j=N p-1 x2 j φ 2 j ( 9
)
The mean error (in the sense of the least squares) between two successive order developments is given by:

2 p ≈ 1 N N k=1 (X p (ξ k ) -X p-1 (ξ k )) 2 ≈ 1 N N k=1 ( Np j=N p-1 +1 xj φ j (ξ k )) 2 (10)
where N is the number of samples.

From the orthogonality property of the gPC (Eq. ( 3)), the last error (Eq. ( 10)) can be written in the following form:

2 p ≈ Np j=N p-1 +1 x2 j φ 2 j = e σ 2 ,p (11)
The equation (Eq. ( 11)) shows a very interesting result; the variance error and the mean least squares error are approximately equal between two successive order developments. It is the basic for finding a criterion for choosing the optimal order of truncation.

The criterion is based on the calculation of the decay rate of the relative error of the variance between two successive order developments and it is given by [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF]:

• γ • γ c 1 c 2 k 2 , k N L 2 k 1 , k N L 1 x 2 x 1 m NES 1 NES 2 (a) NES 1 c h k h , k N L h h 1 (b) NES 2 c h k h , k N L h h 2 (c)
η = e σ 2 ,p σ 2 p ( 12
)
Note that this criterion (Eq. ( 12)) can be estimated directly from the gPC coefficients, it is not necessary to perform additional simulations with the Deterministic Model or with the gPC to calculate these errors. This characteristic represents a major interest of the gPC method 3 Two degrees-of-freedom selfexcited model coupled to NESs

The system considered in this work is composed by the two degrees-of-freedom (DOF) Hultèn's model [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF][START_REF] Hultén | Brake squeal -a self-exciting mechanism with constant friction[END_REF], which represents the primary system, coupled to two identical NESs with masses m h , damping coefficients c h and cubic stiffnesses k N L h . The NESs are attached on the primary system in an ungrounded configuration as shown in Fig. 1.

The equations described the mechanical system are given by :

d 2 x 1 dt 2 + η 1 ω 1 dx 1 dt + ω 2 1 x 1 -γω 2 2 x 2 + ϕ 1 x 3 1 -γϕ 2 x 3 2 + µω 1 dx 1 dt - dh 1 dt + ξ h (x 1 -h 1 ) + ϕ h (x 1 -h 1 ) 3 = 0 d 2 h 1 dt 2 + µω 1 dh 1 dt - dx 1 dt + ξ h (h 1 -x 1 ) + ϕ h (h 1 -x 1 ) 3 = 0 d 2 x 2 dt 2 + η 2 ω 2 dx 2 dt + ω 2 2 x 2 + γω 2 1 x 1 + γϕ 1 x 3 1 + ϕ 2 x 3 2 + µω 1 dx 2 dt - dh 2 dt + ξ h (x 2 -h 2 ) + ϕ h (x 2 -h 2 ) 3 = 0 d 2 h 2 dt 2 + µω 1 dh 2 dt - dx 2 dt + ξ h (h 2 -x 2 ) + ϕ h (h 2 -x 2 ) 3 = 0, (13a) (13b) (13c) (13d)
where h 1 (t) and h 2 (t) (respectively x 1 (t) and x 2 (t)) represent the NESs displacements (respectively the displacements of the primary system),

η i = c i / √ mk i , ω i = k i /m, ϕ i = k N L i /m (with i = 1, 2), = m h /m assuming 0.01 < < 0.1 , ξ h = k h /m, µ = c h / √ mk 1 and ϕ h = k N L h /m.

Problem statement

We focus our analysis on the capacity of the NESs to suppress or mitigate LCOs caused by frictioninduced instabilities. An example of the direct numerical integration of the system with and without NESs is plotted in Fig. 2 which shows the displacements x 1 (t) as a function of the time. Four different values of the friction coefficient are used: γ = 0.16, γ = 0.18, γ = 0.2 and γ = 0.22, the other parameters are set to

ω 1 = 2π100, ω 2 = 2π85, η 1 = 0.02, η 2 = 0.06, ϕ 1 = 10 5 , ϕ 2 = 0, = 0.05, ξ h = 0.001, µ = 0.02, ϕ h = 1.4 • 10 5 . ( 14
)
Moreover, the calculation was performed on 4 seconds in order to be sure that the steady-state has been reached.

Depending on the value of the uncertain parameter, Fig. 2 schows four main types of steady-state regimes generated when two NESs is attached on the primary system: complete suppression of the instability(Fig. 2 We define the amplitude A wN ES 1 of the variables x 1 of the coupled system ( 13) and within a steady-state regime as:

A wN ES 1 = max x SSR 1 (t) -min x SSR 1 (t) 2 , (15)
where x SSR 1 (t) is the times series of the variables x 1 obtained from the numerical integration of the coupled system (13) within the steady-state regime. . This discontinuity corresponds, when γ increases, to the transition from SMR to no suppression regime and separates mitigated regimes and unmitigated regimes.The value of γ at the jump is called mitigation limit (with respect to γ) and denoted γ ml .

The amplitudes A wN ES

In the next section, we present the different methods used in this paper to detect the mitigation limit. The first method considers directly the mitigation as the QoI. The second method considers the amplitudes A wN ES 1 as the QoI, the mitigation is deduced from it.

5 Jump detection methods

The QoI is the mitigation limit

In the first case, the mitigation limit is the quantity of interest. The method is presented in Fig. 4. It involves evaluating the mitigation limit with respect to a deterministic parameter (here we assume that β is determinist) for a given values of friction coefficicient γ. Comparison between time serie x 1 (t) resulting from the numerical integration of the braking system with and without NESs. The set of parameters ( 14) is used.

In order to apply the gPC, the Numerical Experimental Design of the friction coefficient is built using the LHS method (see Sect.2). Then the A wN ES 1 is calculated as a function of the deterministic damping ratio and in a given value of the uncertain friction coefficient. Fig. 5 shows two examples of the variation of A wN ES

The QoI is the amplitude of displacement

In this section, we propose another technique to locate the discontinuity. The amplitude of displacement is now considered as the quantity of interest. Fig. 6 shows the different steps of this method. The damping ratio and the friction coefficient are considered both as uncertain parameters. The Numerical Experimental Design is then built using the LHS method (see Sect.2) with Q simulations. After the construction of the gPC coefficient, the amplitude A wN ES 1 is evaluated using the gPC theory. Then, for a given value of β, the derivative of A wN ES 1 is calculated (see Fig. 7). Finally, the mitigation limit is detected. 14) is used.

Application and results

In this section, the method presented in Sect. 5 is applied to locate the discontinuity in the steady-state amplitude profiles A wN ES 1 (see Eq. ( 15) and Fig. 3) of the system under study [START_REF] Tuan | Design criteria for optimally tuned nonlinear energy sinks-part 1: transient regime[END_REF].

Then the PHSSR of the oscillation by the NESs is estimated. We compute a set of S total samples of the uncertain parameters within the considered uncertain space and following their distribution law. Then, the PHSSR is defined as follows PHSSR = S HSSR S total × 100 [START_REF] George S Fishman | Monte carlo: concepts, algorithms, and applications[END_REF] where S HSSR is the number of samples within the region of uncertain parameters space in which the LCO of the system is mitigated or the system is stable.

Reference study

The reference is evaluated with the two algorithms (see Fig. 4 and Fig. 6). In the case that the mitigation limit is the quantity of interest, 10 000 random iterations of γ and 250 deterministic samples of β are used. The number of simulations is equal to 250 000 (250 x 10 000). In the case that A wN ES 1 is the quantity of interest, 10 000 random iterations of the two uncertain parameters γ and β are used (see Fig. 8). The number of simulations is equal to 10 000 (100 x 100). The result of the reference is shown in Fig. 9. The value of the PHSSR is equal to 92.11% as shown in Tab. 1. Figure 6: Algorithm applied to detect the jump in the case that the quantity of interest is the amplitude of displacement. 6.2 Results of the gPC approach -QoI is the mitigation limit

In this section, the algorithm presented in Fig. 4 is applied to estimate the QoI (mitigation limit) and the PHSSR.

The gPC order is changed from p = 1 to p = 12. As described in Sect. 2.1, the Latin Hypercube Samples (LHS) method is used to build the gPC coefficients using Q simulations where Q = kN p with k = 2.

Fig. 10 shows the results of the proposed algorithm to estimate the mitigation limit. The gPC order p is respectively equal to 2 in Fig. 10(a) and 8 in Fig. 10(b).

The different criteria illustrated in Sect. 2.2 and the comparison between the reference and the proposed method are plotted in Fig. 11. The error curves between the gPC and the reference and also between two successive chaos developments allow to verify that the criterion between two successive orders allows to choose the optimal order of truncation without knowledge of the reference. Moreover, based on the different errors plotted in Fig. 11(a), (b) and (c), p = 6 can be chosen as the optimal order to correctly model the quantity of interest. But for the determination of the PHSSR, the study of the error plotted in Fig. 11(d) shows that the order p = 1 is sufficient. (see ).

Tab. 2 shows the results for p = 1, 2, 6, 8 and 12. In the case where p = 1, the PHSSR is equal to 92% which leads to a relative error of the PHSSR compared with the reference equals to 0.119%. The mean relative error between the reference and the gPC is equal to 3.92% and the number of simulations is equal to 1000. In the case where p = 6, the PHSSR is equal to 92.12% which leads to a relative error of the PHSSR compared with the reference equals to 0.01%. The number of simulations increases from 1000 for p = 1 to 6500 for p = 8. In order to reduce the prohibitive cost of the first method, the algorithm presented in Fig. 6 is now applied to estimate the QoI (the A wN ES

1

) and the PHSSR. The gPC order is changed from p = 8 to p = 20 and Q = 3N p . Fig. 12 shows the results of the proposed algorithm to estimate the mitigation limit. The gPC order p is respectively equal to 8 in Fig. 12(a) and 15 in Fig. 12(b).

The comparison between the reference and the proposed method is plotted in Fig. 13. Fig. 13(a) shows that the mean relative error of the QoI between two successive orders is not correlated with the same error between the reference and the gPC. On the other hand, Fig. 13(c) shows that the relative error of the PHSSR between two successive orders and between the reference and the gPC are correlated. Thus this last criterion and the decay rate allow to choose the optimal order of the gPC, let p = 20.

According to Tab. 3, the PHSSR for p = 20 is equal to 91.66% which leads to a relative error of the PHSSR compared with the reference equals to 0.488%. The number of simulations for the optimal order p = 20 represents 7 times less than the number of simulation needed by the method using the algorithm in Fig. 4 with the optimal order p = 1 and with a relative error of the PHSSR compared with the reference lower than 0.5%. Relative error of the PHSSR (%)

Conclusion

In this paper, a robust modelling of NES used for the brake squeal reduction is presented. For that, a simple two DOF analytical model with uncertain parameters is considered. Related to the dispersion of these uncertain parameters, the coupling of the NES with the model can generate a discontinuity in the evolution of the LCO. Two methods based on the generalized Polynomial Chaos are addressed in order to detect the jump. The first is to consider that the quantity of interest which is a response from the DM as the LCO amplitude and the second is to consider the mitigation limit is the QoI. The aim is to predict the discontinuity and not to obtain an accurate representation of the LCO amplitude. The results of these methods are noted in the following. In the case that the QoI is the LCO amplitude, a polynomial chaos expansion with an order p = 20 must be chosen. In the case that the QoI is the mitigation limit, a polynomial chaos expansion with an order p = 1 must be chosen. Note that the relative error of PHSSR between the reference and the gPC in both these two optimal configurations is less than 0.5%. However, the comparison of the number of simulations allows us to choose which technique is better in terms of computational cost and accuracy: the method that the QoI is the LCO amplitude is the most adequate.

Figure 1 :

 1 Figure 1: (a) Mechanical system with NESs. (b) Zoom on the NES 1 . (c) Zoom on the NES 2 .

  (a)), mitigation through Periodic Response (PR, Fig.2(b)), mitigation through Strongly Modulated Response (SMR, Fig.2(c)) or no mitigation (Fig.2(d)). Bergeot et al.[START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] have been shown these four regimes. Hereafter, mitigated regimes referred to complete suppression, PR and SMR.

1 and A woN ES 1 (

 11 the amplitude of the system without NESs) are plotted as a function of the friction coefficient γ in Fig.3for the set of parameters.The figure highlight a jump (or discontinuity) in the amplitude profile A wN ES 1

  (a) Complete suppression, γ = 0.16 (b) Mitigation: PR, γ = 0.18 (c) Mitigation: SMR, γ = 0.2 (d) No suppression, γ = 0.22

Figure 2 :

 2 Figure 2: Comparison between time serie x 1 (t) resulting from the numerical integration of the braking system with and without NESs. The set of parameters (14) is used.

1 and A woN ES 1

 11 as a function of 250 deterministic samples of β ∈ [1.2 23] (Fig. 5(a) for γ = 0.3 and Fig. 5(b) for γ = 0.4). Then the stability analysis is performed in order to detect the bifurcation point. We denote by β woN ES b the bifurcation point without NESs and β wN ES b with NESs. This study allows us to delimit the stable and the unstable zones of the system. Now the gPC theory is applied, the amplitude A wN ES 1 is evaluated only in the unstable zone when β ∈ [1.2 β wN ES b ]. The mitigation limit is then detected.

Figure 3 :

 3 Figure 3: Amplitudes A wN ES 1

Figure 4 :

 4 Figure4: Algorithm applied to detect the jump when the quantity of interest is the mitigation limit.

Figure 5 : 2 )

 52 Figure 5: The amplitudes of displacement with and without NESs as a function of β: (a) γ = 0.3; (b) γ = 0.4.

Figure 7 :

 7 Figure 7: Example of the variation of the amplitude of displacement with NESs as a function of γ: (a) β = 4.4222; (b) β = 6.7444.

Figure 8 :Figure 9 :

 89 Figure 8: 10 000 deterministic simulations of the amplitude of displacement with NESs as a function of γ and β.

Figure 10 :4Figure 11 :

 1011 Figure 10: The mitigation limit (γ ml ,β ml ) : (a) p = 2; (b) p = 8.

Figure 12 :

 12 Figure 12: The mitigation limit (γ ml ,β ml ) : (a) p = 8; (b) p = 15.

Figure 13 :

 13 Figure 13: Different errors between the reference and the gPC: (a) Relative error of the mean of the QoI; (b) Decay rate %. (c) Relative error of the PHSSR.

Table 1 :

 1 Reference results

		Number of
	PHSSR(%)	simulations using
		algorithm

Table 2 :

 2 Results using the algorithm in Fig.4.

		Number of	Mean			Reltive error
	Order p	simulations from p = 1	relative error reference-	Decay rate(%)	PHSSR(%)	of PHSSR reference-
		to p	gPC(%)			gPC(%)
	1	1000	3.92	-	92	0.119
	2	2000	2.673	0.034	92.06	0.054
	6	6500	0.899	0.002	92.12	0.01
	8	8750	0.788	0.0013	92.13	0.021
	12	12750	0.71	0.001	92.14	0.032
	6.3 Results of the gPC approach -QoI is the amplitude of dis-
	placement

Table 3 :

 3 Results using the algorithm in Fig.6.

					Reltive error
	Order p	Number of simulations	Decay rate(%)	PHSSR(%)	of PHSSR reference-
					gPC(%)
	8	135	2.597	95.78	3.984
	12	321	0.925	94.14	2.203
	15	501	0.498	93.06	1.031
	18	732	1.128	90.53	1.715
	20	899	0.331	91.66	0.488


D T Kieu
email: duc.kieu@insa-cvl.fr

B Bergeot

S Berger

Sparse polynomial chaos expansion for stability analysis of a clutch system with uncertain parameters

Keywords: Stability, Vibration, Clutch, Friction system, Sparse polynomial chaos, Regression methods

In the transmission systems of vehicles, unforced vibrations can be observed during the sliding phase of clutch engagement. These vibrations are due to frictional forces and may generate noise. Several studies have shown that the stability of such friction systems is highly sensitive to parameters (e.g. friction law, damping) which lead to signicant dispersion. Therefore, uncertain parameters must be considered in the stability analysis of a clutch system. In several studies of the literature, the usual generalized polynomial chaos (PC) expansion has been used to study the stability of a clutch system using nonintrusive techniques. However, non-intrusive techniques require a number of model evaluations (i.e. the computational cost) which can become prohibitive when the studied system has a large number of uncertain parameters. To remedy this problem, in this work, we use the sparse polynomial chaos expansion which has been recently developed in reliability domain. The method is compared to the reference Monte Carlo method and with the usual PC expansion in the context of the stability analysis of a clutch system. The results show that the use of the sparse PC allows a remarkable reduction of the computational cost by ensuring a high accuracy compared with the usual PC expansion.

Introduction

Dry friction systems can develop dynamic instabilities related to the friction. In particular, in vehicles with manual transmission systems, these instabilities can be the cause of unwanted vibrations during the sliding phase of the clutch engagement. Several studies have been focused on the mechanisms responsible for these friction-induced vibrations [START_REF] Hervé | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF]. The numerous mechanisms explaining the friction-induced vibration phenomenon are classified into two main families which are related to the tribological aspects of friction systems and to the geometrical and structural properties. The first family explains the instabilities by the variation of the friction coefficient with respect to the relative speed or by a higher static friction coefficient than the dynamic one. The stick-slip is a well-known phenomenon in this context [START_REF] Van De Velde | A new approach of stick-slip based on quasi-harmonic tangential oscillations[END_REF]. The second family attributes the appearance of instabilities to the sprag-slip mechanism and more generally to the mode-coupling phenomenon. In this case, self-excited oscillations may occur even with a constant friction coefficient. Because of the speed range of the measured vibrations, high frequency of the self-excited phenomenon in clutch system cannot be related to stick-slip. Consequently, mode coupling instabilities due to the intrinsic structure of the system are more likely to be responsible for this phenomenon [START_REF] Wickramarachi | Analysis of friction-induced vibration leading to "eek" noise in a dry friction clutch[END_REF]. It is therefore necessary to compute the system eigenvalues which are complex functions because of the friction [4]. The analysis of the real parts sign of the eigenvalues allow to conclude about the stability of the system while the imaginary parts give the frequency of the corresponding mode.

It is known that the dynamic behavior of clutch systems is highly sensitive to design parameters, in particular to the friction coefficient and damping [START_REF] Hervé | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF][START_REF] Centea | Non-linear multi-body dynamic analysis for the study of clutch torsional vibrations (judder)[END_REF]. In addition, the design process can lead to the dispersion of system parameters. Therefore, it is necessary to take into account uncertainties in the system parameters to ensure the robustness of the analysis of friction systems as clutch system and consequently the robustness of the design of this class of systems. The classic Monte Carlo approach which was used to reach this aim requires prohibitive computational cost.

As an alternative, the generalized Polynomial Chaos (gPC) and Multi-Element generalized Polynomial Chaos (ME-gPC) has been used in the past to take into account the uncertainties of the friction coefficient in the study of the dynamic behaviour of friction systems. With gPC [START_REF] Nechak | A polynomial chaos approach to the robust analysis of the dynamic behaviour of friction systems[END_REF][START_REF] Sarrouy | Piecewise polynomial chaos expansion with an application to brake squeal of a linear brake system[END_REF], when the number of uncertain parameters and the order of gPC increase, the number of simulation needed to estimate the PC coefficients becomes extremely large leading to a prohibitive computational cost. To solve this problem, the ME-gPC may be used. In particular, the ME-gPC had been applicated by Trinh et al. [8] for the stability analysis of a clutch system, the conclusion is that the use of the ME-gPC allows a reduction of the computation cost. However, the method is effective with a small number of uncertain parameters (up to 5) but for a higher number of uncertain parameters the computation cost become prohibitive again. To remedy this problem, one considers sparse Polynomial Chaos (sparse PC) expansion. This method is filled out by an iterative procedure allowing to build iteratively a sparse PC expansion while mastering the approximation error [9]. Sparse PC is very efficient for reducing computational cost in reliability analysis of static mechanical structures. This justifies the present work in which sparse PC strategies are used in the context stability analysis of a clutch system model.

The main aim of this study is to investigate the capacity of the sparse PC with Isotropic hyperbolic index sets to study the stability of a clutch system with an increasing number of uncertain parameters compared with the full PC expansion. The main task is to find a harmony between high accuracy and reasonable computational cost.

The paper is organized as follows. In section 2, the stability analysis of dynamical systems with the indirect Lyapounov method is recalled. The Polynomial Chaos theory is briefly presented in section 3 including the description of the gPC, ME-gPC and sparse PC. The section 4 presents the dynamical model of the clutch system under consideration in this study. Finally, the section 5 is dedicated to the comparison of the three polynomial chaos expansions for the stability analysis of the clutch system and to comment on the results.

Stability analysis of dynamic system with the indirect Lyapounov approach

Let us consider the motion equation of a nonlinear dynamic system

ẋ = f (x, u), (1)
where x is the state vector of the system, f is a nonlinear function of x, and u is the uncertain parameter vector.

The Hartman-Grobman theorem (e.g., Wiggins [START_REF] Wiggins | Introduction to applied nonlinear dynamical systems and chaos[END_REF], Chap. 3) states that in the vicinity of a hyperbolic equilibrium point1 , a nonlinear system has the same qualitative stability as does the corresponding linear system. Let the origin x 0 be an equilibrium point for the system (1). The Jacobian matrix of the system (1) is given by [START_REF] Wiggins | Introduction to applied nonlinear dynamical systems and chaos[END_REF]

H(u) = ∂f (x, u) ∂x      x=x 0 . ( 2
)
According to the Lyapunov's indirect method, the stability of x 0 is analyzed by evaluating the eigenvalues λ i (i = 1, ..., n) of the matrix H:

1. The origin x 0 is asymptotically stable if Re(λ i < 0), ∀i ∈ (1, ..., n).

The origin

x 0 is unstable if Re(λ i > 0), ∃i ∈ (1, ..., n).
3 Generalized polynomial chaos and sparse polynomial chaos expansions

Generalized polynomial chaos expansion

The generalized polynomial chaos (gPC) has been proposed by Xiu and Karniakakis [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF] as a generalization of the original Wiener-Chaos expansion. The gPC establishes a separation between the stochastic components of a random function and its deterministic components. From the Wiener theory and the generalizes Cameron-Martin theorem, any second order random process X(ξ) can be expanded in a convergent (in the mean square sense) polynomial function series as

X(ξ) ≈ α∈N r xα φ α (ξ), (3)
where φ α (ξ) are orthogonal polynomial functions, ξ(ξ 1 , ..., ξ r ) is a vector of r independent random variables. The continuous distributions of random variables correspond to the families of orthogonal polynomials (see e.g. Table 1) [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF]. Finally, xα are the polynomial chaos coefficients. In the following, X(ξ) will be called the quantity of interest.

The subscript α is a multi-index, α = {α 1 , ..., α r } in N r . The length of a multi-index α is defined by

|α| = α 1 = r i=1 α i . (4)
The corresponding index set is

(ξ) Support Gaussian Gamma Beta Uniform Hermite Laguerre Jacobi Legendre (-∞,∞) [0,∞) [-1,1] [-1,1] A r,p = α ∈ N r : α 1 ≤ p , (5)
with p the degree of the PC expansion.

The number of terms N p in the gPC expansion, is shown to be dependent on the higher degree p of the polynomials φ α (ξ) and on the stochastic dimension r, i.e. the number of uncertain parameters [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF]

N p = card(A r,p ) = (p + r)! p!r! . (6)
In the regression approach, the coefficients are determined through the minimization of the following least square criterion [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF]

ε 2 reg = Q q=1 X ξ (q) - α∈A r,p xα φ α ξ (q) 2 , (7)
where

ξ (q) = ξ (q) 1 , ..., ξ (q) r
, with q = 1, ..., Q, is the Numerical Experimental Design (NED) and X ξ (q) is the corresponding model evaluations vector. The NED ξ (q) with q = 1, ..., Q may be built using Latin Hypercube Samples method (LHS) [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code[END_REF].

The coefficients are computed as follows

x = φ T (ξ (q) )φ(ξ (q) ) -1 φ T (ξ (q) )X(ξ (q) ), (8)
with φ(ξ (q) ) the matrix defined by

φ(ξ (q) ) =     φ 0 (ξ (1) ) . . . φ Np-1 (ξ (1) ) . . . . . . . . . φ 0 (ξ (Q) ) . . . φ Np-1 (ξ (Q) )     . (9)

Multi-element generalized polynomial chaos expansion

Let ξ( ξ1 , ..., ξr ) denote a ramdom input vector where ξi is independent uniform random variable within the orthogonal interval [-1, 1] r . Next, the space of the random input is decompose into m non-intersecting elements [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF]. A local variables in each element

ζ k (ζ k 1 , ..., ζ k r ) are mapped to a new independent random uniform variables ξk ( ξk 1 , ..., ξk r ) in [-1, 1] r ζ k i = (b k i + a k i )/2 + ξk i (b k i -a k i )/2 i = 1, ..., r; k = 1, ..., m (10)
where a k i , b k i are the lower and upper bounds of the local variables

ζ k i .
Then gPC can be used locally of the k th element, the second random process corresponding to the k th element is given by

X k ( ξk ) ≈ α∈N r xk α φ α ( ξk ). (11)

Sparse polynomial chaos expansion

Isotropic hyperbolic index sets

One defines a hyperbolic truncation from m-norms,

0 < m < 1 [9] A r,p m = α ∈ N r : α m ≤ p , (12)
where

α m = r i=1 α m i 1/m . ( 13
)
Whatever the choice of the value of m-norms, the sequence of nested sets A r,p m (p ∈ N) always converge to the sets N r . Therefore, one defines the isotropic hyperbolic polynomial chaos expansions with the index sets A r,p m

X A r,p m (ξ) = α∈A r,p m xα φ α (ξ). (14)

Error estimates of the polynomial chaos approximations

The strategy consists here in performing an incremental search of the significant terms. For that, it is necessary to define the error estimates of polynomial chaos approximations.

Empirical error. Consider a NED ξ (q) = ξ

(q) 1 , ..., ξ (q) r
, with q = 1, ..., Q, and X ξ (q) is the corresponding model evaluation. X A (ξ) is computed by CP approximation

X A (ξ) ≈ α∈A xα φ α (ξ) (15)
where index set A is a finite non empty subset of N r .

The generalization error is defined by mean of the squared difference between X(ξ

(q) i ) and X A (ξ (q) i ) Err = E X(ξ (q) i -X A (ξ (q) i )) 2 . ( 16
)
In practice, the generalization error may be estimated by the empirical error defined by

Err emp = 1 Q Q q=1 X(ξ (q) ) -X A (ξ (q) ) 2 . ( 17
)
Of common use is the related coefficient of determination R 2 which reads

R 2 = 1 - Err emp V[X] , (18)
where

V[X] is the variance of X ξ (q) V[X] = 1 Q -1 Q q=1 X(ξ (q) ) -X 2 , X = 1 Q Q q=1 X(ξ (q) ). ( 19
)
The empirical error underestimates the generalization error because of the overfitting phenomenon. To avoid this phenomenon, the leave-one-out error, which is a predicted residual sum of squares, is defined below [9].

Leave-one-out error. Let X (-i) A ξ (q) be the chaos polynomial expansion that has been built from the NED ξ (1) , ..., ξ (Q) \ξ (i) , i.e. when removing the i-th observation ξ (i) from the training set ξ (1) , ..., ξ (Q) . The predicted residual is defined as the difference between the model evaluation at ξ (i) and its prediction based on X

(-i) A ξ (q) [9] ∆ (i) = X(ξ (q) ) -X (-i) A (ξ (q) ). ( 20
)
The leave-one-out error so-called the prectited residual sum of squares is defined by

Err LOO = 1 Q Q i=1 ∆ (i) 2 . ( 21
)
It is possible to calculate analytically each predicted residual as follows [9]

∆ (i) = X(ξ (q) ) -X A (ξ (q) ) 1 -h i , (22)
here h i is the i-th diagonal term of the matrix φ(ξ (q) )(φ T (ξ (q) )φ(ξ (q) )) -1 φ T (ξ (q) ). Therefore the leaveone-out error was given by

Err LOO = 1 Q Q i=1 X(ξ (q) i ) -X A (ξ (q) i ) 1 -h i 2 . ( 23
)
The equivalent determination coefficient of R 2 is denoted by S 2

S 2 = 1 - Err LOO V[X] . ( 24
)
The coefficients R 2 and S 2 which are presented in this section are used in an algorithm in the following section to choose an optimal sparse polynomial chaos expansion.

Algorithm applied with the sparse PC expansion

An iterative procedure is presented to build the sparse PC approximation. The objective is to find the best PC coefficients [9].

The algorithm is summarized in 5 basic steps:

Step 1

• Select a numerical experiment design (ξ (q) ), e.g. a random design based on Latin Hypercube Sampling [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code[END_REF]. The model evaluations at the design points are gathered in vector X(ξ

(q) i ).
• Select the values of the algorithm parameters, i.e. the target accuracy S 2 target , the maximal PC degree p max , m-norms of truncation strategy and the cut-off values 1 , 2 .

Step 2

• Initialize the algorithm: p = 0: the truncation index set A 0 = {0} with {0} is the null element of N r .

Step 3: Training step -Enrichment of the pc basis (valid for any degree p ∈ [1, ..., p max ])

⇒ Forward step (Addition step):

-Gather the candidate terms in a set C r,p m = α ∈ N r : p -1 ≤ α m ≤ p , add each candidate term to set A r,p-1 m one-by-one and compute the CP expansion coefficients by regression using Eq. 8 and the associated determination coefficient R 2 in each case.

-Retain eventually those candidate terms that lead to a significant increase in the coefficient R 2 , i.e. which significantly decrease the empirical error Err emp , and discard the other candidate terms. Let A r,p+ m be the final truncation set at this stage.

⇒ Backward step (Elimination step):

-Remove in turn each term in A r,p+ m of degree strictly less than p. In each case, compute the PC expansion coefficients and the associated determination coefficient R 2 .

-Eventually discard from A r,p+ m those terms that lead to an insignificant decrease in R 2 , i.e. a negligible increase of the empirical error Err emp . Let A r,p m be the final truncation set.

Step 4

• Verify the conditioning of the regression information matrix. If it is poor, i.e. the size of the NED (ξ (q) ) is smaller than 2card(A p m ) with isotropic index sets, a enrichment of the numerical experiment design is done. In this case, the truncation set A is reset to {0} and the basis enrichment procedure is restarted using nested Latin Hypercube designs [START_REF] Wang | Adaptive Response Surface Method Using Inherited Latin Hypercube Design Points[END_REF].

Step 5: Test step • Stop if either the leave-one-out error S 2 is less than the target S 2 target or the order of the PC basis is equal to p max .

The algorithm applied with the isotropic index sets for building up a sparse polynomial chaos expansion is presented in figure 1.

Analytical model of the clutch system

The lumped parameter model of figure 2 that used in this study was defined by P. Wickramarachi [START_REF] Wickramarachi | Analysis of friction-induced vibration leading to "eek" noise in a dry friction clutch[END_REF] considering in addition the internal damping of the clutch disk and the pressure plate. This model with 6 degrees of freedom has been chosen because it is sufficient and efficient to study the instabilities generated by mode couplings and has been validated by experimental results. In this model, the contact

M • Ü + C • U + K • U = 0 (25)
with

U = θ x θ y Z A Z B Z C Z D T , (26)
M = diag I x I y Mp 4 Mp 4 Mp 4 Mp 4 , (27)
K =          r 2 (k A + k B + 4k f ) -µlr (k C + k D ) r (k A + 2k f ) -r (k B + 2k f ) µlk C -µlk D -µlr (k A + k B ) r 2 (k C + k D + 4k f ) -µlk A µlk B r (k C + 2k f ) -r (k D + 2k f ) r(k A + 2k f ) 0 k A + 2k f 0 -k f -k f -r(k B + 2k f ) 0 0 k B + 2k f -k f -k f 0 r(k C + 2k f ) -k f -k f k C + 2k f 0 0 -r(k D + 2k f ) -k f -k f 0 k D + 2k f          , (28)
C =          r 2 (c A + c B + 4c f ) -µlr (c C + c D ) r (c A + 2c f ) -r (c B + 2c f ) µlc C -µlc D -µlr (c A + c B ) r 2 (c C + c D + 4c f ) -µlc A µlc B r (c C + 2c f ) -r (c D + 2c f ) r(c A + 2c f ) 0 c A + 2c f 0 -c f -c f -r(c B + 2c f ) 0 0 c B + 2c f -c f -c f 0 r(c C + 2c f ) -c f -c f c C + 2c f 0 0 -r(c D + 2c f ) -c f -c f 0 c D + 2c f          , (29)
k A = γ 1 k p /4; k B = k p /( 4γ 1 )
;

k C = γ 2 k p /4; k D = k p /(4γ 2 ), (30)
N A = k A (Z A -rθ y ) N B = k B (Z B + rθ y ) (31)
N C = k C (Z C + rθ x ) N D = k D (Z D -rθ x ), F A = µN A F B = µN B (32)
F C = µN C F D = µN D ,
where r = (r 1 + r 2 )/2 with r 1 and r 2 are the minimum and maximum sliding radi, µ is the friction coefficient which is assumed constant because of the the high slip speed (<700 rev/min) and l is the thickness of the pressure plate. The nominal values of the parameters are: k p = 16 MN/m; k f = 7 MN/m; γ 1 = 0.9; γ 2 = 0.8; r 1 = 75 mm; r 2 = 120 mm; l = 12.5 mm;

c A = c B = c C = c D = 4N m -1 s -1 ; c f = 0.1N m -1 s -1 [3].

Application and results

The objective of this section is to study the stability of a clutch system with uncertain parameters. These studies are done in the lumped parameter presented in section 4. There are eight parameters which can be uncertain: µ, k p , k f , γ 1 , γ 2 , r 1 , r 2 and l. The uncertain parameters are considered as independent random and uniform in the intervals [0.95*V n , 1.05*V n ] with V n the corresponding nominal values; µ is independent random and uniform in the interval [0, 0.5]. The nominal values of the parameters have been given in section 4.

Reference study: Monte Carlo on the initial model

In this section, the stability of the clutch system is investigated using the Monte Carlo on the initial model (MCIM). 10,000 eigenvalues λ i (i = 1, ..., n) of the Jacobian matrix (Eq. ( 2)) of the clutch system are calculated directly from a set of 10,000 independent random samples.

In the eigenvalues of the clutch system, there are four non-zeros eigenvalues (modes 1, 2, 3 and 4 with the complex conjugates). Modes 3 et 4 are always decoupled and stable and the coalescence phenomenon can occur between the modes 1 and 2. Therefore, the stability of the clutch system only depends on these two modes.

Figure 3 shows the real and imaginary parts of mode 1 (blue) and mode 2 (red) and highlights the coalescence phenomenon of eigenvalues described hereafter. From µ = 0, the imaginary parts (frequencies) of these two modes are separated and tend to come closer when the coefficient of friction increases. As for the real parts, they are both negative, i.e. the system is stable. Then, the imaginary parts of the two modes overlap, it is the so-called mode coalescence phenomenon. Beyond this point, imaginary parts are almost equal. From the Hopf bifurcation point, if µ increases, the real part of one mode is positive and the system becomes therefore unstable.

Sparse PC method for the stability analysis of a clutch system

In this section, the distribution of the eigenvalues of the clutch system is rebuild by coupling the MC method with the polynomial chaos expansion methods presented in section 3.

Table 2 presents the comparison between the results obtained with gPC, ME-gPC, Sparse PC and MCIM. The results using gPC and MC-gPC were computed by Trinh et al. [8]. The results show that gPC can be used in the stability analysis of the clutch system with a small number of uncertain parameters r (up to 5). With a high value of r, the number of Direct Computations (DC) of the initial system needed to determine the PC coefficients became too high and the relative errors of the proportion stability between gPC and MCIM are higher than 20% from r = 6. Te use of the ME-gPC is effective for r ≤ 7. However, with r = 5 → 7, the numbers of DC of the complete system are rather high (> 3000), that means the computational costs are expensive. With r = 8, the number of DC to determine the coefficients of the PC expansion exceed the limited number (10,000) for both methods.

With Sparse PC, authors compute five sets of NED to analyze the influence of the choice of the NED on the Sparse PC method. A target accuracy S 2 target = 0.999 and a maximal PC degree p max = 15 are used. In table 2, the results of each of the five NED are not all depicted, only the minimal and maximal values are shown in square brackets. The results show that the use of this method reduce remarkably the number of DC compare with gPC and ME-gPC, specially for r ≥ 5. Moreover, the method ensures a high accuracy compared with gPC. The relative errors of the proportion of stability using Sparse PC are lower than 5% for r = 3 → 5 than 7% for r = 6 → 8

For r = 8, gPC and ME-gPC are not effective because the number of DC exceeds the limited number whereas the number of DC using Sparse PC lower than 3000. It means a saving of 70% of the computational cost with a high accuracy (< 7%).

Figure 4 shows the eigenvalue real parts of modes 1 and 2 using MCIM (blue) and using Sparse PC (red) (for r = 3, 4, 5, 6, 7, 8). The results computed with Sparse PC are close to the results obtained via MCIM. The Hopf bifurcation points with Sparse and MCIM are also close. Therefore, Sparse PC ensures a high accuracy of the results.

Conclusion

The computational cost of gPC and ME-gPC in studying the stability analysis of a clutch system can become prohibitive with a large number of uncertain parameters (r ≥ 5). While using Sparse PC with Isotropic hyperbolic index sets allows a remarkable reduction of the computational cost by ensuring a high accuracy compared with these common PC expansion. With number of uncertain parameters r = 8, gPC and ME-gPC were not feasible because of the limited number of DC. Whereas using Sparse PC isotropic hyperbolic index sets reduce even more than 70% of the computational cost with a high accuracy (the relative errors of the proportion stability are lower than 7% with Sparse PC anisotropic hyperbolic index sets.

In conclusion, Sparse polynomial chaos is efficient for high number of uncertain parameters (r = 3 → 8 in this study). It reduces remarkably the computational cost by ensuring a high accuracy compared with the common polynomial chaos expansion.

Figure 1 :

 1 Figure 1: The algorithm applied with the Isotropic hyperbolic index sets for building up a sparse polynomial chaos expansion

Figure 2 :

 2 Figure 2: Analytical model of the clutch system

Figure 3 :

 3 Figure 3: (a) Real part and (b) imaginary parts of the eigenvalues λ 1 and λ 2 for the modes 1 and 2

Figure 4 :

 4 Figure 4: Real parts of mode 1 and 2 using MCIM and Sparse PC

Table 1 :

 1 Correspondence between the families of orthogonal polynomials and random variable Random variable ξ The polynomial family φ α

If the eigenvalues of the Jacobian matrix of the nonlinear system evaluated at the fixed point have non zero real parts the fixed point is hyperbolic.