Stability analysis of a clutch system with multi-element generalized polynomial chaos

– 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 (friction law, damping) which lead to signiﬁcant dispersions. Therefore, uncertain parameters must be considered in the stability analysis of a clutch system. This paper investigates the ability of the multi-element generalized polynomial chaos to take an increasing number of uncertain parameters into account in the stability analysis of a clutch system: it focuses on accuracy, on the criterion for the choice of the order of truncation, on the criterion for the choice of the number of elements and on computation time. The objective is to propose a low-cost model of high accuracy, compared with the Monte Carlo method.


Introduction
In vehicles with manual transmission systems, unforced vibrations can be observed during the sliding phase of clutch engagement.Several studies have focused on the mechanisms responsible for these self-excited frictioninduced vibrations [1].Various mechanisms have been identified to explain the friction-induced vibration phenomenon.They are classified into two main families which are related to the tribological aspects of friction systems and to the geometrical and structural properties.Low frequency phenomena such as judder (10−20 Hz) can often be attributed to tribological properties [2].However, high frequency phenomena cannot be related to stick-slip behaviour because of the speed range of the measured vibrations.Consequently, mode coupling instabilities due to the intrinsic structure of the system are more likely to be responsible for this phenomenon [3].It is thus necessary to compute the system eigenvalues which are complex functions because of the friction [4].The real parts allow the study of the system stability and the imaginary parts give the frequency of the corresponding mode.
Numerous studies have also demonstrated that the dynamic behaviour of dry friction systems is highly sensitive to design parameters, in particular to the friction coefficient and damping [5][6][7].For instance, Chevillot studied the effects of damping and in particular the destabilization paradox in an aircraft braking system [8].Hervé studied the effects of friction and damping on the stability of clutch systems; his study focused on the destabilization paradox, but the model was a lumped model with only two degrees of freedom (DOF) [1,9].The friction and damping effects on the mode coupling phenomenon in the finite element (FE) model of a vehicle brake squeal were shown in the studies of Fritz [10,11].In these studies, in order to reduce the calculation time which may be prohibitive because of the high number of DOF, damping is taken into account in a modal approach.It is therefore necessary to study the influence of the dry friction coefficient and damping and, in particular, the destabilization paradox of the instability of the clutch system in a modal approach.
Moreover, the manufacturing process has been shown to lead to dispersions in the design parameters.That is why the dispersion of the uncertain parameters must be taken into account to ensure the robustness of the analysis of friction systems and thus the robustness of the design of this class of systems.The Monte-Carlo (MC) approach which is classically used to reach for this purpose requires prohibitive calculation times, especially for a system with a high number of DOF.As an alternative, the polynomial chaos formalism has been proposed to take account of the uncertainties of the friction coefficient in the study of the dynamic behaviour of friction systems [12,13].But these studies were carried out on the model of a braking system with two DOF and there were only one or two uncertain parameters.In addition, the limitation of generalized polynomial chaos (gPC) requires a large amount of calculations when the number of stochastic modes is large.To solve this problem, a method based on multi-element generalized polynomial chaos (ME-gPC) is proposed in several papers.Nevertheless, these studies again just consider one or two uncertain parameters [4,[12][13][14][15] or are theoretical and only apply to a mathematical system of differential equations [4,16,17].
The main objective of this paper is therefore to investigate the abilities and the limitations, in terms of accuracy and computational costs, of gPC and ME-gPC to take account of an increasing number of uncertain parameters (up to 8) in modelling the eigenvalues of a clutch system and in analyzing the system stability.The challenge is to find a compromise between accuracy and computation time.The aim is now to assess the maximum number of uncertain parameters which can be considered with the gPC expansions, in particular in the present study of the stability of a clutch model.The latter involves a physical model which has been validated experimentally.The paper focuses on the criteria for choosing the order of truncation and the number of elements, independently of the results of the MC method.The results of the methods based-on the gPC expansions are compared with those obtained with the classic MC approach for validation.The final objective is to propose a low-cost method of high accuracy, in order to avoid the computational difficulties of the classic MC method.
This paper is organized as follows: Section 2 presents the methods for the analysis of the system stability and Section 3 is devoted to gPC and ME-gPC formalisms; the friction system is described in Section 4; the deterministic study of the influence of parameters on the mode coupling of a clutch system in a modal approach is presented in Section 5, with a focus on the destabilization paradox; the results of the robust stability analysis are given in Section 6, followed by a conclusion in Section 7.

Analysis of the stability of a dynamic system 2.1 Stability analysis of the eigenvalues
Consider the motion equation of a nonlinear dynamic system: where X represents the instantaneous state of the system (its coordinates in the phase space), the upper dot denotes the time derivation, f is a function of X and d stands for uncertain parameters.According to the Hartman-Grobman theorem, the linearization of Equation ( 1) in the vicinity of the equilibrium state X e (d 0 ) preserves its non-marginal stability nature.Therefore, the determination of the stability nature of equilibriums only requires the knowledge of the linearized motion equation in their vicinity in most cases.Because of the form of the solutions, the stability nature of X e (d 0 ) is expressed by the eigenvalues λ(d 0 ) of the Jacobian Df (X e (d 0 ), d 0 ).So, in the classic MC procedure, to analyze the stability of a system which has uncertain parameters , samples are first generated following the probabilistic support of those parameters, then the eigenvalues λ(d) corresponding to each sample are calculated.This sampling-based method is known to be time-consuming since a large number of samples are needed to ensure reasonable accuracy with high confidence.The resulting computing costs are exorbitant since the eigenvalues of the system must be calculated for each sample, which is difficult, especially for systems with numerous DOF.Therefore, the gPC formalism can be used instead of the classic MC procedure.

Modal approach
The real parts of complex eigenvalues help to study the system stability and the imaginary parts give the frequency of the corresponding mode.The modal approach presented here allows the calculation of the eigenvalues of the coupling modes.Using the modal approach, the number of retained modes can be truncated, which reduces the number of DOF and computation time.It also helps to better understand the influence of damping with its management by mode.
Consider the linearized equation of motion of a dynamic system: where M, C, K are respectively the mass, damping and stiffness matrices.u is the displacement vector and the dot denotes the derivative with respect to time.Because of the friction, the stiffness matrix has special properties: with K s the structural stiffness matrix, K f the asymmetric stiffness matrix caused by the friction, and μ the friction coefficient.Equation ( 2) defines an eigenvalue problem in which the eigenvalues λ and eigenvectors ψ must be determined.To solve Equation (2), two approaches are possible: the direct approach or the modal approach.Here, the latter will be detailed, which involves projecting Equation (2) on the real modal basis of the system before extracting the complex eigenvalues.This requires the calculation of the real eigenvalues considering only the symmetric part of the stiffness matrix, in order to find the real frequencies and the modal basis of the frictionless system associated with Equation (4).
After the extraction of n mode modes of Equation ( 4), let Ω = diag(ω 1 . . .ω nmode ) be the diagonal matrix gathering the corresponding frequencies, and γ the vector of the displacement coordinates in the undamped and frictionless modal basis.The displacement vector in Equation ( 2) can be expressed as a linear combination of the undamped and frictionless modes.Once transformed in the frequency domain and expanded using the orthogonality relations, Equation (2) can be rewritten as: where D and Λ f are the projections of C and K f on the undamped and frictionless modal basis.I is the identity matrix and s the Laplace parameter.
For each value of μ, Equation (5) defines a problem with complex eigenvalues of dimension (n mode × n mode ) which is much smaller than the first dimension (2n DOF × 2n DOF ) of Equation ( 2).This size is small enough to allow the problem to be treated using numerical computation software.The modal approach is much more effective than the direct approach in terms of computation time.

Base theories of generalized polynomial chaos
Generalized polynomial chaos (gPC) establishes a separation between the stochastic components of a random function and its deterministic components.In the dynamic system (1), if the uncertain parameters are uniform, all the eigenvalues λ i (i = 1, ..., n) are also random functions.According to the Askey scheme [2], the Legendre polynomials are best suited to deal with uniform uncertainties, so the random eigenvalues are given by: φ j (ξ) are the Legendre polynomials [2].ξ (ξ 1 , ..., ξ r ) is a vector of r independent random variables, ξ i is distributed uniformly within the orthogonality interval [-1,1].
The truncation order N P is shown to be dependent on the polynomial chaos order p and stochastic dimension r denoting the number of uncertain parameters: The computation of λ i is then turned into the problem of finding the coefficients λi,j of its truncated expansion.These coefficients, called stochastic modes, can be computed with the non-intrusive spectral projection (NISP) or the regression methods.The principal advantage of these methods lies in the fact that no modification is performed on the system; only the calculation of the eigenvalues of the completed clutch system for a limited number of samples is required.In the regression method, the stochastic modes λi,j are computed by minimizing the least square criterion: where ξ (q) (q = 1, ..., Q) are well known Gauss collocation points [2].In practice, the Q points are selected from a set of roots of the Legendre polynomials [2,17].If the same order is chosen for all random variables, Q is calculated as follows: More details on the techniques used to estimate the stochastic modes can be found in [12].

Statistical characteristics of the gPC
The main advantage of the expansion of gPC is the direct estimation of the statistical characteristics of stochastic processes when stochastic methods are estimated.Indeed, the moments of the first and second orders are given by: ⎧ ⎪ ⎨ ⎪ ⎩ λi = λi,0 where λi and σ 2 i are the mean value and the variance of the ith eigenvalue [4].

Criteria for the choice of the optimal truncation order
The higher is the truncation order, the better is the accuracy of gPC, but also the higher are the amount of stochastic modes and calculation time.It is therefore necessary to determine the optimal truncation order.This section presents the criteria for the choice of the optimal truncation order in terms of accuracy and computation time.
Let λ i,p ,λ i,p−1 be the eigenvalues vectors which are successively calculated with the two gPC developments of orders p, p − 1.
In the literature, the optimal truncation order is often chosen according to the mean error ε 2 i,p in the L 2 -sense between two developments of successive orders [2,12,14]: with N the number of samples.
It must be noted that it is necessary to calculate the eigenvalues obtained with the gPC in order to estimate the mean error ε 2 i,p in the L 2 -sense.The error of the variance between two developments of successive orders is: From the orthogonality of gPC, we have shown (Appendix 1): From Equation ( 13) it can be seen that the mean error, which is found equal to the error of variance, can be directly computed from the stochastic modes instead of performing simulations with the initial system (Eq.( 5)) (MC) or gPC expansion.This characteristic property is the basis for finding a criterion for the selection of the optimal truncation order.
The criterion considers the decay rate of the relative error between two developments of successive orders [4].This decay rate is defined by: The first criterion for choosing the optimal truncation order is defined as: where θ is called the threshold related to the choice of the truncation order.
In addition, as discussed in Section 2.1, the use of gPC is really effective if the number of direct calculations (DC) of the complete system for the construction of stochastic modes is lower than the number of calculated samples with the MC method.
Therefore, the second criterion for choosing the optimal truncation order concerns the number Q of DC for the construction of the stochastic modes: where N is the number of calculated samples with the MC method to ensure reasonable accuracy with great confidence in the calculation of the stability proportion for a given random space [18].The new feature and the advantage of these criteria to estimate the optimal order is that they are based on the knowledge of stochastic modes and it is not necessary to estimate the eigenvalues in a large number of samples with gPC, contrary to other criteria in the literature [2,12,14].

Decomposition of the random space
As mentioned in the previous section, the approximation error with gPC may be controlled by increasing the value of the order p.However, the size of the stochastic system, represented by N p , and the number of DC to construct the stochastic modes, represented by Q, will increase rapidly.gPC will therefore be difficult to implement in a complex system with numerous DOF and a large number of uncertain parameters r.
This has led to the design of the layout of multielement generalized polynomial chaos (ME-gPC).This technique is mainly based on the decomposition of the random space into m non-intersecting elements [4,16].
In practice, the local variables in each element ζ k are expressed in terms of independent random uniform variables ξk in [-1,1] r .
where a k i , b k i are the limits of the uniform intervals of the kth element.
The GPC expansion of the random process corresponding to the kth element and to the ith eigenvalue of the system is given by:

Statistical characteristics of the ME-gPC
The approximate local mean value and variance of the ith eigenvalue in the kth element of gPC for a truncation order p are: The approximate global mean value and variance for a truncation order p can then be written respectively as: The approximations of the stochastic processes are given by: where J k denotes the probability measurement [4,16].

Criteria for choosing the optimal number of elements
As mentioned in Section 3.1.3,the limitation of gPC is driven by the need to restrict the number of DC for the construction of the stochastic modes, thus the polynomial order cannot be increased significantly in an arbitrary way.However, ME-gPC can overcome this difficulty.In this approach, the error contribution of the kth element can be decreased by a factor J k (dictated by the size of the element), even if the approximated local error is high.In ME-gPC, a small order of gPC can be used for each element since the degree of local disturbance has been reduced.
The selection criterion for the optimal number of elements also considers the decay rate of the relative error in each element [4].This decay rate is defined by: For the decomposition of the random space, two factors are considered: the decay rate η k and the factor J k .A random element will be divided into two equal parts when the following condition is satisfied: When the random elements become smaller (i.e.J k becomes smaller), the value of η k satisfying the criterion will be allowed to be greater.Thus, the criterion relaxes the restriction on the accuracy of the local variance for small elements as the error contribution of small random elements will be dictated by their size.
The second criterion for the choice of the optimal number of elements is based on the determination of the most sensitive parameter.The sensitivity of each parameter is defined by: where the index i, p indicates the mode whose p order associated polynomial function of gPC only depends on parameter ξ i .All the parameters that satisfy Equation (25) will be divided into two random elements equally in the next time step, while all other parameters will remain unchanged: The third criterion for choosing the optimal number of elements is the number of DC needed for the construction of the stochastic modes.Equation (16) shows that the number of DC in an element is (p + 1) r , so the total number of DC for the construction of the stochastic modes in all the m elements is: This paper proposes a new strategy based on the decomposition algorithm for each element according to the above criteria.This strategy is presented in Figure 1.
According to the decomposition algorithm, each element which does not satisfy criteria ( 23) and (25) will be divided into two equal parts, so the maximum number of elements after the ith division level is 2 ni and the element size in the ith division level is calculated with: where J k is the probability measurement.Thus, the maximum number Q of DC needed for the construction of the set of stochastic modes must satisfy: Minimum element size J k min is therefore defined by: A random element will be divided into two equal parts when conditions ( 23) and (25) and the following condition are satisfied:

Squeal model of the clutch system
There are not industrial models in the literature for stability analysis of clutch systems.Wickramarachi has proposed a squeal model involving six DOF which is sufficient and efficient to study the mode coupling instability and which has been validated with experiments [3].In the model, the contact between the friction disc (1) and the flywheel ( 2) is created at points A , B , C , D by a progressive spring k p which is split into four springs k A , k B , k C and k D (Fig. 2).In order to consider the nonlinear characteristic of the progressive spring, the stiffnesses k A and k B are respectively divided and multiplied  The linear model of the undamped clutch system can be expressed as follows: see equation (34) above

Study of the influence of parameters on the mode coupling of a clutch system in a deterministic approachdestabilization paradox
The aim of this section is to study the influence of the most significant parameters on the stability of the clutch system.The studies are carried out in the modal approach with a focus on the destabilization paradox of damping.

Influence of friction coefficient μ
In the clutch system, there are four non-zero eigenvalues (with their conjugates) which depend on the values μ, k p , k f , (modes 1, 2, 3 and 4).Modes 3, 4 are always decoupled and stable; the coalescence phenomenon occurs between modes 1 and 2. The system stability only depends on these two modes.
Figure 3 shows the real parts and imaginary parts of mode 1 (red) and mode 2 (blue).The curves show the well-known phenomenon of mode coupling or coalescence between mode 1 and mode 2. From μ = 0, the imaginary parts (frequencies) of the two modes are separated and they tend to come closer when the friction coefficient increases.When the real parts of the modes are negative, the system is stable.Then, the two modes reach the same imaginary part at a point called the point of coalescence.Beyond this point, the imaginary parts of the two modes are equal, but their real parts become non-zero and opposite.One mode becomes unstable because its real part becomes positive from the Hopf bifurcation point and the other mode remains stable because its real part remains negative.So, the system becomes increasingly unstable if μ increases.

Influence of damping on the stability of a clutch system
The aim of this section is to calculate the eigenvalues with the modal approach in order to study the influence of modal damping on the stability of a clutch system.The modal damping matrix, of size n modes × n modes , is of the form D = diag(0, . . .0, d 1 , d 2 , 0,. . ., 0).D represents the damping projection on the modal basis of the frictionless system.d 1 and d 2 are respectively the damping coefficients for modes 1 and 2. Table 1 shows 16 cases of study divided in four major groups: undamped, equally damped, with slight and large differences in damping.It can be seen that if d i increases, the real parts of the eigenvalues decrease, while the imaginary parts are not influenced.The real parts cross the axis 0 for a higher value of the friction coefficient.Therefore, the higher the damping of the clutch is, the more stable it becomes.This conclusion is confirmed in Table 1 with an increase in the Hopf bifurcation point μ 0 (cases 2 to 4).

Non-equally damped coalescence
Non-equally damped coalescence has been studied in the following cases: slight differences in damping values, large differences in high and small damping values (Table 1).The results are compared to the case of equally damped coalescence with small damping.First, we consider the cases with slight differences (cases 5 to 6) and the cases with large differences in high damping values (cases 11 to 16).Such damping distributions deeply alter the coalescence patterns, as displayed in Figures 5 and 6.As expected, an increased damping mainly tends to shift the curves towards the negative values of the real parts.The clutch behaviour is thus improved in terms of stability, as illustrated in Table 1 by an increase in the Hopf  bifurcation point μ 0 .Nevertheless, the coupling patterns are far more complicated than before.
First of all, this damping distribution induces a gap between the real parts of the two modes at μ = 0.Then, the main change occurs in the vicinity of the coalescence point, where the transition is less sharp than previously.As explained by Hoffmann [19], a "smoothing effect" of the curves with respect to the friction coefficient is observed, both for real parts and imaginary parts.As a consequence, the real parts seem to start to split at a lower friction coefficient value than in the undamped situation.Moreover, the higher the damping difference is, the larger the gap between the imaginary parts of the two modes is.Even if the coalescence patterns are more complicated than previously, they feature a noteworthy advantage.Indeed, since the imaginary and real parts are different, for each friction value, the unstable mode can now be clearly identified and tracked with respect to the friction coefficient value.So far, an increase in damping seemed to stabilize the system.It shifted the curves towards the negative real part values and so, induced an increase in the Hopf bifurcation point μ 0 .However, this is not always the case.Indeed, a difference in damping has two main effects on the coalescence curves: a "lowering effect" for the real part values, and a "smoothing effect" in the vicinity of the coalescence point.The first effect was predominant in the curves shown previously.But with a large difference in damping values between the two modes (cases 7 to 10, Table 1), the second effect may prevail, as illustrated in Figure 7.The two damped curves plotted in Figure 7 respectively present a ratio of 2.5 and 5 between the damping values of the two modes.In these cases, the smoothing effect induces a significant drop in μ 0 , which leads to a drop in stability, as shown in Table 1.Thus, the addition of a certain amount of damping to the system, unevenly distributed on the two modes, leads to instability.This phenomenon is called the destabilization paradox of damping.

Stability maps
This section aims at synthetizing the effects of damping on the clutch stability.It can be inferred from the   previous section in which it has been shown that the ratio between the damping values of the two modes highly influences the clutch stability.
In this analysis, the damping of mode 1 is set to a nonzero constant value, while the damping of mode 2 is chosen so that the ratio d 2 /d 1 ranges from 0 to 5. Figures 8a  and 8b display the stability maps (μ, d 2 /d 1 ) respectively with d 1 = 4 and d 1 = 10.For each couple of values, the clutch stability is assessed through the sign of the real parts of the eigenvalues.The red area is unstable while the blue area is stable.The abscissas of the border between the two areas are by definition the Hopf bifurcation points μ 0 .When the two damping values are different and of low level (d 1 = 4), the stability becomes significantly worse (Fig. 8a).The optimal case in terms of stability appears to be the equally damped case (d 2 /d 1 = 1).When the two damping values are different and of high level (d 1 = 10), the stability area increases when damping d 2 grows for d 2 /d 1 < 1 and also d 2 /d 1 > 1 (Fig. 8b).
Thus, the destabilization paradox of damping can be generalized to a number of values of μ when damping is unevenly distributed with either small or high damping values.This section therefore shows that the dynamic behaviour of the clutch system is highly nonlinear and sensitive to design parameters.Moreover, it highlights and confirms the destabilization paradox in a modal approach for a clutch system with 6 DOF.This paradox has been studied by Hervé et al. on clutch model with only two DOF [1].

gPC and ME-gPC approach to the stability analysis of a clutch system with six DOF
The objective of this section is to investigate the capabilities of the polynomial chaos approach, in terms of accuracy and computational costs, to model the eigenvalues of the system and to analyze the system stability with a high number of uncertain parameters.The results of the gPC method are compared with those obtained with the classic MC approach in order to validate the gPC method.

Uncertain parameters
In the clutch system with six DOF shown in Section 4, there are generally eight uncertain parameters: μ, k p , k f , γ 1 , γ 2 , r 1 , r 2 and l.In this section, the uncertain parameters k p , k f , γ 1 , γ 2 , r 1 , r 2 and l are considered as random and uniform in the intervals [0.95V n , 1.05V n ] with V n the corresponding nominal values; μ is random and uniform in the interval [0, 0.5].The nominal values of the parameters have been given in Section 4.
In order to analyze the stability of the clutch system with the gPC approach, the advantages of the modal approach (shown in Sect.5) are used to calculate the eigenvalues.Here, the study focuses on modes 1 and 2 (and their conjugates).The modal damping values are equal (d 1 = d 2 = 6).
In this section, new studies aiming at highlighting the efficiencies of the gPC and the ME-gPC for large numbers of uncertain parameters will be presented.The selection criteria of the optimal truncation order and of the optimal number of elements will be illustrated.

gPC approach
The stochastic modes are calculated with gPC in eight cases which correspond to the number of uncertain parameters r = 1 → 8. First, the eigenvalues are modelled with gPC to determine the stochastic modes with the regression method using the Q Gauss collocation points and the modal approach.Then, the optimal truncation order is chosen with the validation of the relevance of the first criterion (Eq.( 15)).All the results are compared to those obtained with the MC method with the complete system.
Figure 9 shows the evolutions of the real and imaginary parts of mode 1 (circle line) and mode 2 (cross line) with DC and gPC with r = 1 and order p = 19.The evolutions are nearly similar.
Figure 10 shows the evolutions of the mean errors (in the L 2 -sense) of the real parts of the eigenvalues, between DC and gPC (red line), between two successive orders using directly the real parts of gPC (Eq.( 11)) (blue line) and between two successive orders gPC using the stochastic modes of gPC (Eq.( 12)) (black line).The blue and black lines are very close, which confirms that Equation ( 13) is verified, that is to say, the mean error (in the L 2 -sense) is equal to the variance error between two developments in successive orders of gPC.   Figure 11 shows the evolution of the decay rate η i (Eq.( 14)) (blue line) and the evolution of the relative variance error of the real parts of mode 1 between DC and gPC (red line).Figure 12 shows the evolutions of the relative errors of the stability proportion respectively between DC and gPC (red line), and between two successive orders using gPC (blue line).The curves show the same evolution as the "decay rate" curve and the errors are reduced in the same way if p increases; therefore, the choice of the "decay rate" criterion is also found relevant for the evaluation of the quality of gPC expansion in the estimation of the system stability.The optimal order p is selected by increasing the value of p until the percentages of the relative errors are lower than the thresholds (Eq.( 15)), the second criterion (Eq.( 16)) being verified.
In this study, the threshold value has been set to θ = 0.02 which corresponds to the relative variance error  of 0.04%.The models based on gPC were constructed with r uncertain parameters (r = 1 → 8) and the proportion of stability was calculated with these models and one variation of μ.The stability of the equilibrium point is analyzed for each of the N samples generated according to the probabilistic law considered for uncertain parameters (N = 10 000 in this study).N was set to ensure a confidence level of 99% with an accuracy margin of 1.25% for the stability proportion [18].The values of the optimal orders obtained for the studies are shown in Table 2.
Figure 13 shows the evolution of the real parts of mode 1 and mode 2 using gPC (blue) and the DC (red) with r = 3.The Hopf bifurcation points obtained with the two methods are close.The errors are acceptable in these cases (see Table 1).Figure 14 shows the results using gPC with r = 7.In this case, the errors are quite significant, showing that the optimal order with the fixed threshold cannot be found for the first criterion (here θ = 0.02).Indeed, the second criterion (Eq.( 16)) avoids increasing the order (here N = 10 000).
Table 2 shows the comparison of the results of the stability analysis between the use of gPC and the initial system respectively.The strategy for finding the optimal order consists in stopping the calculations when one of the two criteria is met -either the decay rate η i < 4e − 4  or the number of calculations ≥ 10 000.The optimal order can be found up to r = 4.For a greater value of r (r > 5), no order p can be found to satisfy the first criterion (Eq.( 17)), p is therefore selected according to the second criterion (Eq.( 16)).In Table 2, the relative errors of the results of the stability analysis are lower than 9% up to r = 5, and higher than 20% from r = 6.Thus, the use of these two criteria gives a good estimate of the stability proportion up to r = 5, but is not sufficient for r > 5.In conclusion, these results show that the application of gPC in the stability analysis of the clutch system with five uncertain parameters is effective.The choice of the optimal truncation order with the "decay rate" criterion which is independent of the results obtained with the MC method is appropriate.Furthermore, in the case of a small r (r < 4), the number of DC of the complete system (for the construction of the stochastic modes) is markedly lower than the number of calculations with the MC method.Therefore, the computational costs are greatly reduced.

ME-gPC approach
Section 6.2 shows the effectiveness of gPC for the analysis of the stability of a clutch system with six DOF in cases with a small number of uncertain parameters r (up to r = 5).With a high r value, the number of DC of the complete system is too high and gPC is no more effective.That is why a ME-gPC approach is proposed to overcome such issues.
An optimal number of elements m is usually selected by increasing the value of m while the relative errors (Eq.( 22)) in all the elements verify Equation (23).However, this method is not effective if ME-gPC is used with regular intervals.Indeed, when the criteria are not verified, all the elements are divided; the stochastic modes must then be reconstructed for all the elements, while this is not necessary for the elements for which the chaos expansion is already efficient.Indeed, the quality of the results obtained with ME-gPC may be highly variable according to the intervals.For example, in Figure 13 (in the case of gPG with r = 3), the results calculated with gPC are very different from the results calculated with DC for the interval μ ∈ [0.075, 0.125] and very close in some intervals (e.g.μ ∈ [0, 0.075] or μ ∈ [0.125, 0.5]).gPC is already effective in these intervals.That is why a better algorithm is necessary for finding the optimal m.This algorithm consists in the decomposition method of the random space, which is described with the criteria (Eqs.( 23), ( 25) and (30)) in Section 3.2.3 (Fig. 1).
In the next study, we choose θ 1 equal to 0.02 and 0.01 successively which corresponds to a relative variance error of 0.04% and 0.01% respectively, θ 2 = 0.99 which corresponds to the decomposition of the random space for the most sensitive parameter and N = 10 000 as the maximum number of DC of the complete system for the criterion in Equation (30).The resulting optimal values of the numbers of elements are shown in Table 3.
Figures 15-17 show the evolutions of the real parts of modes 1 and 2 using DC (red), with gPC (blue) and using ME-gPC with the optimal number of elements m (r = 1, 3, 5 with the corresponding p fix ) (green).The results calculated with ME-gPC are closer to the results obtained via DC than the results calculated with gPC.The Hopf bifurcation points with ME-gPC and DC are also close.Therefore, ME-gPC has improved the accuracy of the results.
Table 3 presents the comparison of the results successively obtained with DC, with gPC and the optimal   order p, and with ME-gPC and the optimal number of elements m.For each value of r and p, the calculated results correspond to a given value of θ 1 .The results naturally show that if θ 1 is smaller, the number of elements is higher and the error of the stability proportion is lower.The optimal value for this threshold θ 1 is 0.02.
With ME-gPC, the accuracy of the results has improved significantly.The relative errors of the results of the stability analysis are lower than 3% up to r = 7.Moreover, with a small value of r (up to 3), the number of DC of the complete system (for the construction of stochastic modes) using ME-gPC is higher than the number of calculations using gPC.With r = 4 and 5, the number of DC using ME-gPC is lower than the number of calculations with gPC.With r = 6 and 7, gPC is not effective for the stability analysis (the relative error is higher than 20%) because of the limited number of DC, but ME-gPC is effective in these cases (the relative error is lower than 3%).ME-gPC will therefore be more effective with a medium value of r (r = 4 to 7).With r = 8, because of the limited number of DC, the gPC and ME-gPC methods are not effective for the stability analysis of a clutch system.It is worth noticing that the results in this paper were obtained using a conventional PC equipped with an Intelcore I7-3520M cpu 2.9 GHz.The orders of magnitude of the computational times ranges from 1s up to 10 s for the calculation of the stochastic modes, and are respectively lower than 1s for the calculation of the eigenvalues with the gPC developments and lower than 1s for the calculation of the eigenvalues with the 6 DOF clutch model.

Conclusion
The dynamic behaviour of the studied clutch system is highly nonlinear and highly sensitive to design parameters.This paper confirms the destabilization paradox in a modal approach for a clutch system.In order to avoid the drawbacks of the classic MC method whose cost is prohibitive for industrial systems, gPC and ME-gPC have been proposed.The ability of gPC and ME-gPC to take an increasing number of uncertain parameters into account in the stability analysis of a clutch system has been investigated.Two criteria for the choice of the truncation order have been used; they consist in the relative error of the variance of the eigenvalues between two developments of successive orders and the maximum number of DC.Similarly, three criteria for the choice of the number of elements have also been used; the first one is the relative variance error of the eigenvalues between two developments of successive orders for each element, the second one is a criterion which helps to choose the most sensitive parameters for which the intervals must be divided; the third one is the minimal size of each element, which depends on the maximum number of DC.Polynomial chaos can be efficient for up to seven uncertain parameters and computation time is significantly reduced.A higher number of uncertain parameters -using the necessary truncation orders and numbers of elements to obtain a high accuracy of polynomial chaos -requires a greater number of DC of the complete system than the number needed in the classic MC approach.The results show that the application of polynomial chaos in the stability analysis of

Fig. 1 .
Fig. 1.ME-gPC algorithm used to split the random dimension.

Figure 4
Figure4shows the results corresponding to three low values of equally distributed damping coefficientsd 1 = d 2 .It can be seen that if d i increases, the real parts of the eigenvalues decrease, while the imaginary parts are not influenced.The real parts cross the axis 0 for a higher value of the friction coefficient.Therefore, the higher the

Table 1 .
Damping and Hopf bifurcation point.

Table 2 .
Comparison of results between DC and gPC.

Table 3 .
Comparison of the results respectively obtained using DC, gPC with the optimal order p, ME-gPC with the optimal number of elements m.