Issue 
Mechanics & Industry
Volume 20, Number 1, 2019



Article Number  104  
Number of page(s)  17  
DOI  https://doi.org/10.1051/meca/2019003  
Published online  28 March 2019 
Regular Article
Stability analysis of a clutch system with uncertain parameters using sparse polynomial chaos expansions
INSA CVL, Univ. Orléans, Univ. Tours, LaMé EA 7494, 3 Rue de la Chocolaterie, CS 23410, 41034 Blois Cedex, France
^{*} email: baptiste.bergeot@insacvl.fr
Received:
31
August
2018
Accepted:
3
January
2019
In vehicle transmission systems, frictional forces acting during the sliding phase of the clutch engagement may produce unwanted vibrations. The prediction of the stability of a clutch system remains however a laborious task, as the parameters which have the highest impact on the stability, such as the friction law or the damping, lead to significant dispersions and must be considered as uncertain in such studies. Nonintrusive generalized polynomial chaos (gPC) expansions have already been used in this context. However, the number of deterministic model evaluations (i.e. the computational cost) required to compute the PC coefficients becomes prohibitive for large numbers of uncertain parameters. The sparse polynomial chaos, recently developed by Blatman and Sudret, may overcome this issue. In this paper, the method has been applied to the stability analysis of a clutch system owning up to eight uncertain parameters. Comparisons with the reference Monte Carlo method and classic full PC expansions show that sparse PC expansions allow substantial computational cost reductions while ensuring a high accuracy of the results.
Key words: Stability / vibration / clutch / friction system / sparse polynomial chaos / regression
© D.T. Kieu et al., published by EDP Sciences 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Dynamic instabilities resulting in inconvenient vibrations and noise may appear in dry friction systems. Such phenomena are for instance regular during the sliding phase of the clutch engagement in vehicles equipped with a manual transmission system. Various mechanisms may lie behind those instabilities [1]. In a first category which gathers the mechanisms related to tribological considerations, the instabilities are linked to the dependence of the friction coefficient to the relative speed, or to the difference between the dynamic and static friction coefficients, the latter being usually slightly higher. The stickslip phenomenon is for instance typical of this family [2,3]. Geometrical or structural aspects are implied in the second category of mechanisms, where the instabilities are ascribed to the spragslip or to the modecoupling mechanisms and may occur even when the friction coefficient remains constant.
In a clutch system, the high velocities at stake in the observed phenomenon make it more likely due to a mode coupling instability than to stickslip [4]. The evaluation of the system eigenvalues is hence a necessary step to analyze the stability, which is related to the sign of their real parts (the eigenvalues being complex in the presence of friction [5]). The imaginary parts are linked to the mode frequencies. Past studies have shown that the eigenvalues are extremely sensitive to design parameters, such as the friction coefficient and the damping [1,6–9]. The effect of those two parameters has been investigated in particular in [1,6], in which a lumped model of a clutch system owning two degrees of freedom (DOF) is considered, and where the destabilization paradox is given a special interest. The difficulty in the stability analysis of a clutch system lies in the fact that those design parameters are characterized by a large dispersion, due notably to the high variability of the operating conditions [10], and to the manufacturing process. That dispersion must therefore be taken into account in the analysis of the dynamic behavior. A stochastic approach seems ideal to ensure the robustness of such an analysis.
In the classic Monte Carlo (MC) approach, the number of required deterministic model (DM) evaluations (i.e. the computational cost) may be huge. If the method is efficient for small models, it becomes too expensive for larger systems owning a high number of DOF, such as finite element models used in an industrial context.
More effective methods, such as the generalized polynomial chaos (gPC) and the multielement generalized polynomial chaos (MEgPC), have been developed in this context. They consist in expanding the quantity of interest in series of polynomials of uncertain variables. The coefficients of the PC expansion, which are deterministic and denoted PC coefficients, have to be evaluated, using either intrusive or nonintrusive methods. The intrusive approach requires to modify the DM to determine the PC coefficients. In the nonintrusive strategy, the latter are evaluated from a number of simulations of the deterministic system. The gPC has been used to study the dynamic behavior of friction systems with two DOF and with only one or two uncertain parameters. For instance, Nechak et al. analyzed the stability of a break system respectively with a direct Lyapunov approach coupled with an intrusive gPC [11], and an indirect Lyapunov approach associated with a nonintrusive gPC [12]. Moreover, Sarrouy et al. [13] used an intrusive gPC to study the stability of the same simplified 2 DOF break model with an indirect Lyapunov method. The same method was applied in [14] to analyze the dispersion of the eigenvalues and eigenvectors of a finite element rotor with a limited number of uncertain parameters. The nonintrusive approach, which is more suitable in an industrial context, has been chosen in this paper. It has been applied for instance in [15], along with Wiener–Haar expansions, to analyze the break squeal from a finite element model of a brake system with one uncertain parameter. In [16] and [17], hybrid approaches associating the nonintrusive gPC with respectively the inclusion function based on Chebyshev polynomials and the kriging formalism are considered to predict frictioninduced instabilities with 2 and 3 uncertain parameters defined either by a probability law or by intervals. A Padé expansion based on gPC has been developed in [18] to obtain the dynamic response in the frequency domain of a 2 DOF system with up to 2 uncertain parameters. The drawback of the method lies in the fact that the number of DM evaluations required to compute the PC coefficients increases with the gPC order p and the number of uncertain parameters r. The computational cost may therefore also become excessive for high p and/or r values. In that cases, the MEgPC may be used as an alternative. For instance, the MEgPC appeared efficient to predict the frictioninduced vibrations of a nonlinear uncertain dry friction system in a study of Nechak et al. [19]. More recently, Trinh et al. [20] have applied the MEgPC to analyze the stability of a clutch system, and achieved a substantial reduction of the computational cost in comparison with the standard gPC. However, the method was proved effective for a limited number of uncertain parameters, the computational cost being again excessive when more than 5 uncertain parameters were considered.
To overcome this issue, sparse polynomial chaos (sparse PC) expansions, i.e. which contain a low number of nonzero coefficients compared to the gPC or MEgPC expansions, may be considered. They are inspired by the socalled sparsityofeffects principle [21], which states that most models are principally governed by main effects and loworder interactions, hence limiting the number of useful polynomials coupling different uncertain parameters. In this paper, three different sparse PC expansions are presented. The first one is the sparse PC with lowrank index sets, which has been first defined in [22,23]. This method restricts the number of interacting parameters in the polynomials that couple several parameters. The second and third methods are respectively the sparse PC with isotropic hyperbolic index sets and the sparse PC with anisotropic hyperbolic index sets [24], which limit the orders of the terms in the retained polynomials, keeping high degree polynomials for a single parameter and lower degrees for the polynomials expressing the interactions between several parameters. In the third variant, the use of additional weights favors the parameters that are seen the most sensitive. For each method, an iterative procedure allows to further reduce the retained index set on the basis of error computations. The algorithms are presented in [25,26] for sparse PC with lowrank index sets, and in [24] for sparse PC with hyperbolic index sets.
The studies cited above show that sparse PC expansions are very efficient for reducing the computational cost in reliability analyzes of static mechanical structures. Developing sparse PC strategies in the dynamic field therefore constitutes an original challenge. The aim of this paper is to investigate the capabilities and limitations of the sparse PC expansions to study the stability of a specific dynamic system, namely a clutch system, with an increased number of uncertain parameters compared to that traditionally considered with gPC and MEgPC expansions. The underlying issue is to find a compromise between the desired accuracy level and the corresponding computational cost. The results of these methods are compared with those obtained with the classic MC approach for validation.
The rest of the paper is organized as follows: Section 2 presents the retained lumped model of the clutch system along with the indirect Lyapunov approach to perform the stability analysis. The MC method is described in Section 3. Section 4 exposes the gPC, MEgPC and the strategies related to sparse PC expansions. The application of the proposed methods to analyze the stability of the clutch system and the corresponding results are finally given in Section 5.
2 Analytical model of the clutch system
2.1 Equations of motion
The parameters of the clutch system lumped model considered in this paper are presented in Figure 1. This 6 DOF model, originally defined by Wickramarachi [4] and validated experimentally, is suitable for the study of modecoupling induced instabilities as it involves a sufficient but limited number of DOF. The model has been completed to take into account additional damping linked to the clutch disk and the pressure plate [20]. The details are recalled hereafter for the sake of clarity. As shown in Figure 1, the friction surface of the clutch system is connected to the flywheel at points Aʹ, Bʹ, Cʹ and Dʹ through stiffnesses k _{ A }, k _{ B }, k _{ C } and k _{ D } which are related to the stiffness k _{ p } of the cushion spring. In order to consider the nonlinearity of the clutch friction disk stiffness value with respect to the position of the pressure plate, the stiffnesses k _{ p }/4 at locations Aʹ and Bʹ are respectively multiplied and divided by a factor γ _{1}; a second factor γ _{2} is similarly used at locations Cʹ and D ʹ, as detailed in equation (1) below [4]$${k}_{A}={\gamma}_{1}{k}_{p}/4,\text{}\text{\hspace{1em}}{k}_{B}={k}_{p}/(4{\gamma}_{1})\text{,}$$(1a) $${k}_{C}={\gamma}_{2}{k}_{p}/4,\text{\hspace{1em}}{k}_{D}={k}_{p}/(4{\gamma}_{2})\text{.}$$(1b)
The internal damping of the clutch disk is taken into account through four damping coefficients c _{ A }, c _{ B }, c _{ C } and c _{ D } placed at the same locations than the above stiffnesses. The contact surface between the friction disk and the pressure plate is a ring of internal radius r _{1} and external radius r _{2}. The points Aʹ, B ʹ, C ʹ and D ʹ are supposed located at the average radius $\overline{r}}\mathrm{}=({r}_{1}+{r}_{2})/2$ of that ring. The pressure plate of mass M_{p} and thickness 2l is modeled by four masses M _{ p }/4 connected by a bending stiffness k _{ f } and a damping coefficient c _{ f } (see Fig. 1 on the right). The points A, B, C and D represent the projections of A ʹ, B ʹ, C ʹ and D ʹ on its average surface, whereas the points E, F, G, H are fixed points of the flywheel. The 6 DOF are the internal rotations θ _{ x }, θ _{ y } around the fixed axes x and y which describe rigidbody rotations (wobbling modes) of the pressure plate, and the translation movements Z _{ A }, Z _{ B }, Z _{ C }, Z _{ D } of points A, B, C, D which depict the first two nodaldiameter bending modes along the z axis.
The normal spring forces ${N}_{{A}^{\prime}}$, ${N}_{{B}^{\prime}}$, ${N}_{{C}^{\prime}}$, ${N}_{{D}^{\prime}}$ and the corresponding friction forces ${F}_{{A}^{\prime}}$, ${F}_{{B}^{\prime}}$, ${F}_{{C}^{\prime}}$, ${F}_{{D}^{\prime}}$ at points A ʹ, B ʹ, C ʹ and D ʹ depend on the above DOF as follows:$${N}_{{A}^{\prime}}={k}_{A}({Z}_{A}+{\displaystyle \overline{r}}{\theta}_{x}),\text{\hspace{1em}}{N}_{{B}^{\prime}}={k}_{B}({Z}_{B}{\displaystyle \overline{r}}{\theta}_{x})\text{,}$$(2a) $${N}_{{C}^{\prime}}={k}_{C}({Z}_{C}+{\displaystyle \overline{r}}{\theta}_{y}),\text{\hspace{1em}}{N}_{{D}^{\prime}}={k}_{D}({Z}_{D}{\displaystyle \overline{r}}{\theta}_{y})\text{,}$$(2b)and$${F}_{{A}^{\prime}}=\mu {N}_{{A}^{\prime}},\text{\hspace{1em}}{F}_{{B}^{\prime}}=\mu {N}_{{B}^{\prime}}\text{,}$$(3a) $${F}_{{C}^{\prime}}=\mu {N}_{{C}^{\prime}},\text{\hspace{1em}}{F}_{{D}^{\prime}}=\mu {N}_{{D}^{\prime}}\text{.}$$(3b)
In the above equations, the friction coefficient μ is assumed constant because of the high velocity range (>700 rev/min). The resulting motion equations of the clutch system model are written in matrix form as$$\mathbf{M}\stackrel{}{\mathbf{U}}+\mathbf{C}\dot{\mathbf{U}}+\mathbf{K}\mathbf{U}=\mathbf{0},$$(4)with$${\mathbf{U}=[{\mathit{\theta}}_{\mathit{x}}\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\mathit{\theta}}_{\mathit{y}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\text{\hspace{0.17em}Z}}}_{\mathit{A}}{\mathit{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Z}}}_{\mathit{B}}{\mathit{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Z}}}_{\mathit{C}}{\mathit{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Z}}}_{\mathit{D}}]}^{\mathit{T}}\mathrm{.}$$(5)
In equation (4), the mass matrix is given by$$\text{M}=\text{diag}\left(\left[\begin{array}{cccccc}\hfill {I}_{x}\hfill & \hfill {I}_{y}\hfill & \hfill \frac{{M}_{p}}{4}\hfill & \hfill \frac{{M}_{p}}{4}\hfill & \hfill \frac{{M}_{p}}{4}\hfill & \hfill \frac{{M}_{p}}{4}\hfill \end{array}\right]\right)\text{,}$$(6)where I _{ x }, I _{ y } are the respective moments of inertia around the axes x and y.
The stiffness and damping matrices are defined respectively as $$\text{K}=\left[\begin{array}{cccccc}\hfill {{\displaystyle \stackrel{}{r}}}^{2}\left({k}_{A}+{k}_{B}+4{k}_{f}\right)\hfill & \hfill \mu l{\displaystyle \stackrel{\u203e}{r}}\left({k}_{C}+{k}_{D}\right)\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({k}_{A}+2{k}_{f}\right)\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({k}_{B}+2{k}_{f}\right)\hfill & \hfill \mu l{k}_{C}\hfill & \hfill \mu l{k}_{D}\hfill \\ \hfill \mu l{\displaystyle \stackrel{\u203e}{r}}\left({k}_{A}+{k}_{B}\right)\hfill & \hfill {{\displaystyle \stackrel{\u203e}{r}}}^{2}\left({k}_{C}+{k}_{D}+4{k}_{f}\right)\hfill & \hfill \mu l{k}_{A}\hfill & \hfill \mu l{k}_{B}\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({k}_{C}+2{k}_{f}\right)\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({k}_{D}+2{k}_{f}\right)\hfill \\ \hfill {\displaystyle \stackrel{\u203e}{r}}({k}_{A}+2{k}_{f})\hfill & \hfill 0\hfill & \hfill {k}_{A}+2{k}_{f}\hfill & \hfill 0\hfill & \hfill {k}_{f}\hfill & \hfill {k}_{f}\hfill \\ \hfill {\displaystyle \stackrel{\u203e}{r}}({k}_{B}+2{k}_{f})\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {k}_{B}+2{k}_{f}\hfill & \hfill {k}_{f}\hfill & \hfill {k}_{f}\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}({k}_{C}+2{k}_{f})\hfill & \hfill {k}_{f}\hfill & \hfill {k}_{f}\hfill & \hfill {k}_{C}+2{k}_{f}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}({k}_{D}+2{k}_{f})\hfill & \hfill {k}_{f}\hfill & \hfill {k}_{f}\hfill & \hfill 0\hfill & \hfill {k}_{D}+2{k}_{f}\hfill \end{array}\right]$$(7)
and $$\text{C}=\left[\begin{array}{cccccc}\hfill {{\displaystyle \stackrel{\u203e}{r}}}^{2}\left({c}_{A}+{c}_{B}+4{c}_{f}\right)\hfill & \hfill \mu l{\displaystyle \stackrel{\u203e}{r}}\left({c}_{C}+{c}_{D}\right)\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({c}_{A}+2{c}_{f}\right)\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({c}_{B}+2{c}_{f}\right)\hfill & \hfill \mu l{c}_{C}\hfill & \hfill \mu l{c}_{D}\hfill \\ \hfill \mu l{\displaystyle \stackrel{\u203e}{r}}\left({c}_{A}+{c}_{B}\right)\hfill & \hfill {{\displaystyle \stackrel{\u203e}{r}}}^{2}\left({c}_{C}+{c}_{D}+4{c}_{f}\right)\hfill & \hfill \mu l{c}_{A}\hfill & \hfill \mu l{c}_{B}\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({c}_{C}+2{c}_{f}\right)\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}\left({c}_{D}+2{c}_{f}\right)\hfill \\ \hfill {\displaystyle \stackrel{\u203e}{r}}({c}_{A}+2{c}_{f})\hfill & \hfill 0\hfill & \hfill {c}_{A}+2{c}_{f}\hfill & \hfill 0\hfill & \hfill {c}_{f}\hfill & \hfill {c}_{f}\hfill \\ \hfill {\displaystyle \stackrel{\u203e}{r}}({c}_{B}+2{c}_{f})\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {c}_{B}+2{c}_{f}\hfill & \hfill {c}_{f}\hfill & \hfill {c}_{f}\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}({c}_{C}+2{c}_{f})\hfill & \hfill {c}_{f}\hfill & \hfill {c}_{f}\hfill & \hfill {c}_{C}+2{c}_{f}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \stackrel{\u203e}{r}}({c}_{D}+2{c}_{f})\hfill & \hfill {c}_{f}\hfill & \hfill {c}_{f}\hfill & \hfill 0\hfill & \hfill {c}_{D}+2{c}_{f}\hfill \end{array}\right]\text{.}$$(8)
The nominal values of the parameters, taken from [4], are: k _{ p } = 16 × 10^{6} kg·s^{−2}; k _{ f } = 7 × 10^{6} kg·s^{−2}; γ _{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 } = 4 kg·s^{−1}; c _{ f } = 0.1 kg·s^{−1}.
Fig. 1 Analytical model of the clutch system. On the left two rotations θ _{ x }, θ _{ y } (blue arrows) around the x and y axes, on the right four translation movements Z _{ A }, Z _{ B }, Z _{ C }, Z _{ D } (green arrows) of points A, B, C, D in the z direction. 
2.2 Stability of the trivial equilibrium solution
The stability of the trivial equilibrium solution U _{0} = 0 is investigated by calculating the eigenvalues λ of the equations of motion (4) solving the following characteristic equation$$\mathrm{det}\left({\lambda}^{2}\text{M}+\lambda \text{C}+\text{K}\right)=0\text{.}$$(9)
Because the system (4) is a 6 DOF secondorder system it has 12 eigenvalues. Indeed, the characteristic equation (9) is a 12th degree polynomial equation with respect to λ.
The Lyapunov's indirect method relates the stability of ${\text{U}}_{0}$ to the signs of the eigenvalues λ _{ i } (i = 1,…, 12):

${\text{U}}_{0}$ is asymptotically stable if ∀i ∈ (1,…, 12), Re (λ _{ i }) < 0;

${\text{U}}_{0}$ is unstable if $\exists \mathit{i}\text{\hspace{0.17em}}\in \text{\hspace{0.17em}}(1,\mathrm{\dots},12)/\mathit{R}\mathit{e}\mathit{(}{\mathrm{\lambda}}_{\mathit{i}}\mathit{)}>0.$
The squeal stability analysis described above is a complex eigenvalue analysis (CEA). It is known that CEA may lead to an underestimation or an overestimation of the number of unstable modes due to the limitations of the linear assumption [27]. Indeed, a complete nonlinear analysis shows that additional unstable modes can appear during transient vibrations. However, the determination of the linear stability of an equilibrium point using CEA is the first step in the global analysis (i.e. squeal starts under linear conditions around an initial equilibrium point). To be sure that all fundamental frequencies and associated harmonics or harmonic combinations eventually contribute to friction induced vibrations, it would be necessary to calculate the selfexcited vibrations via an integration scheme. It is however not the purpose of the present paper, in which the goal is only to determine if the system is stable or unstable.
In the stochastic study presented in the following sections, the eigenvalues λ _{ i } (i = 1,…, 12), solutions of the characteristic equation (9), are the quantities of interest (QoI).
3 Monte Carlo method
In the classic MC procedure, to analyze the stability of a system with uncertain parameters, some samples are generated following the probabilistic support of these parameters. The corresponding eigenvalues are then computed for each sample. This method is hereafter referred to as the MC method applied to the DM (MCDM). Besides, in the MC method coupled with the PC expansion, the eigenvalues λ _{ i } are computed with the PC expansion for each sample. In the MC simulation, if N independent random samples are used, the means of the QoI (here the eigenvalues) converge to the rate of O (N ^{−1/2}), whatever the dimension of problem [28]. In this paper, 10 000 independent random samples are generated following the probabilistic support of the parameters to analyze the stability of the clutch system. The polynomial chaos theory is briefly reviewed in the following section.
4 Generalized polynomial chaos and sparse polynomial chaos
4.1 Generalized polynomial chaos
The gPC, originally proposed by Xiu and Karniadakis [29], consists in expanding any second order process X (ξ) depending on r independent random variables (ξ _{1},…, ξ _{ r }) = ξ into the series$$X(\xi )={\displaystyle \sum _{\alpha \in {\mathbb{N}}^{\mathrm{r}}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}(\xi )\text{,}$$(10) where ϕ _{ α } (ξ) are orthogonal polynomials which represent the stochastic components of the process, and ${{\displaystyle \overline{x}}}_{\alpha}$ are the PC coefficients that account for the deterministic components of the process. The vector ξ ∈ [−1, 1] ^{ r } is linked to the vector of the physical uncertain parameters $\zeta =({\zeta}_{1},\mathrm{\dots},{\zeta}_{r})\in {\prod}_{i=1}^{r}[{a}_{i},{b}_{i}]$ through the following relation$${\zeta}_{i}=\frac{{b}_{i}+{a}_{i}}{2}+{\xi}_{i}\frac{{b}_{i}{a}_{i}}{2},\text{\hspace{1em}}i=1,\dots ,r\text{,}$$(11)where a _{ i } and b _{ i } denote the limits of the uniform interval corresponding to the ith parameter.
The Wiener theory as well as the generalized Cameron–Martin theorem [30] state that the series is convergent in the mean square sense. According to the Askey scheme, if the distributions of the uncertain parameters are uniform (as considered in the following), the polynomial functions ϕ _{ α } are most suitably obtained from Legendre polynomials, as detailed in Table 1 [29,31,32].
In practice, the QoI X(ξ) are approached by a truncated expansion involving a finite number of terms N _{ p }:$$X(\xi )\approx {\displaystyle \sum _{\alpha \in {\mathcal{A}}^{r,p}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}(\xi )\text{,}$$(12)where p is the order of the PC expansion and α ={ α _{1},…, α _{ r } } a multiindex in ${\mathbb{N}}^{\mathit{r}}$, of length$$\left\right\alpha {}_{1}={\displaystyle \sum _{i=1}^{r}}{\alpha}_{i}\text{.}$$(13)
The index set used in the truncated expansion (12) is then defined as$${\mathcal{A}}^{r,p}=\{\alpha \in {\mathbb{N}}^{\mathit{r}}:\left\alpha \right{}_{1}\le p\}\text{.}$$(14)
Computing the QoI X is therefore turned into the problem of finding the coefficients ${{\displaystyle \overline{x}}}_{\alpha}$ of the truncated gPC expansion (Eq. (12)). The number of terms N _{ p } to evaluate is linked to the order p and to the number of uncertain parameters r as [29]$${N}_{p}=\text{card}({\mathcal{A}}^{r,p})=\frac{(p+r)!}{p!r!}\text{.}$$(15)
As mentioned above, these coefficients may be determined either from Galerkin intrusive methods or from nonintrusive methods. In an industrial context, the latter are more appropriate because they do not require any modification of the DM; the PC coefficients are built from a finite number of values of the QoI X, computed from numerical simulations of the DM. The nonintrusive spectral projection (NISP) and the regression methods are the most used nonintrusive methods. The NISP method requires simulations at the Q = (p + 1)^{ r } Gauss collocation points to construct the PC coefficients. With the regression method, a minimum of Q = N _{ p } simulations are necessary, which may be performed in practice at the Q = (p + 1)^{ r } Gauss collocation points or at points chosen with a Latin Hypercube Samples (LHS) method [33]. In the latter case, Q = kN _{ p } simulations (with k a small integer usually equal to 2, 3 or 4) are used.
Within the regression framework, the evaluation of the coefficients results from the minimization of the following criterion [34]$${\epsilon}_{\text{reg}}^{2}=\mathrm{}{\displaystyle \sum _{q=1}^{Q}}{[X\left({\xi}^{\mathrm{}(q)\mathrm{}}\right){\displaystyle \sum _{\alpha \in {\mathcal{A}}^{r\mathrm{},p}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}\left({\xi}^{(\mathrm{}q)\mathrm{}}\right)}^{2},\mathrm{}$$(16)where ${\xi}^{(q)}=({\xi}_{1}^{(q)},\mathrm{\dots},{\xi}_{r}^{(q)})$ (with q = 1,…, Q) denotes the Numerical Experimental Design (NED), that is a set of Q vectors of parameter values generated from the probabilistic support of the parameters; X (ξ ^{(q)}) denotes the vector of the corresponding DM evaluations. The PC coefficients are finally computed as$$\overline{x}}={({\varphi}^{T}({\xi}^{(q)})\varphi ({\xi}^{(q)}))}^{1}{\varphi}^{T}({\xi}^{(q)})X({\xi}^{(q)})\text{,$$(17)with $\mathit{\varphi}({\xi}^{(q)})$ the matrix defined by$$\varphi ({\xi}^{(q)})=\left(\begin{array}{ccc}\hfill {\varphi}_{0}({\xi}^{(1)})\hfill & \hfill \dots \hfill & \hfill {\varphi}_{{N}_{p}1}({\xi}^{(1)})\hfill \\ \hfill \vdots \hfill & \hfill \ddots \hfill & \hfill \vdots \hfill \\ \hfill {\varphi}_{0}({\xi}^{(Q)})\hfill & \hfill \dots \hfill & \hfill {\varphi}_{{N}_{p}1}({\xi}^{(Q)})\hfill \end{array}\right)\text{.}$$(18)
The method to choose the optimal truncation order is described in [20], where a criterion based on the decay rate of the relative error between two expansions of successive orders is used. In practice, the decay rate is defined as the ratio between the error of the variance between two expansions of orders p − 1 and p, and the variance of the expansion of order p.
If the number of uncertain parameters and the order p of the gPC expansion are high, the number of PC coefficients and therefore the necessary number of simulations to build them become quickly prohibitive. Strategies to reduce this simulation number are consequently necessary.
Correspondence between the families of orthogonal polynomials and the distributions of the random variables.
4.2 Multielement generalized polynomial chaos
This approach consists in breaking down the random space into m nonintersecting elements, in the aim of reducing the number Q of simulations needed to evaluate the PC coefficients by using a PC expansion with a low order p in each element [35].
As in the gPC method, a change of variables is performed, in each of the m elements, from the vector of the local uncertain variables ${{\displaystyle \overline{\zeta}}}^{k}\in {\prod}_{i=1}^{r}[{a}_{i}^{k},{b}_{i}^{k}]$, to the vector of independent uniformly distributed random variables ${{\displaystyle \overline{\xi}}}^{k}\in {[1,1]}^{r}$, that is$${{\displaystyle \overline{\zeta}}}_{i}^{k}=\frac{{b}_{i}^{k}+{a}_{i}^{k}}{2}+{{\displaystyle \overline{\xi}}}_{i}^{k}\frac{{b}_{i}^{k}{a}_{i}^{k}}{2},\text{\hspace{0.22em}}\text{\hspace{0.22em}}\text{\hspace{0.22em}}i=1,\dots ,r;\text{\hspace{0.22em}}\text{\hspace{0.22em}}\text{\hspace{0.22em}}k=1,\dots ,m\text{.}$$(19)
In this context, the random process corresponding to the kth element is given by the local gPC expansion$${X}^{k}({{\displaystyle \overline{\xi}}}^{k})\approx {\displaystyle \sum _{\alpha \in {\mathcal{A}}^{r,p}}}{{\displaystyle \overline{x}}}_{\alpha}^{k}{\varphi}_{\alpha}({{\displaystyle \overline{\xi}}}^{k})\text{,}$$(20)while the approximation of the stochastic process is finally expressed using probability measurements J _{ k } [35] as$$X(\xi )\approx {\displaystyle \sum _{k=1}^{m}}{\displaystyle \sum _{\alpha \in {\mathcal{A}}^{r,p}}}{{\displaystyle \overline{x}}}_{\alpha}^{k}{\varphi}_{\alpha}({{\displaystyle \overline{\xi}}}^{k}){J}_{k}\text{.}$$(21)
The detailed procedure to build an MEgPC expansion can be found for instance in [20], where the method was applied to analyze the stability of a clutch system. The results obtained in [20] will be shown for comparison purposes in the last section, but the present paper is dedicated to the use of sparse PC expansions as a strategy to reduce the computational costs.
In the context of MEgPC, a low order of the PC expansion in each element can be kept, as the accuracy of the results is mainly driven by the number m of elements chosen to decompose the random space. Strategies to define an optimal number m are described in [20], where an iterative procedure using three criteria is proposed. For a given p value, the first criterion is again based on the decay rate of the relative error between two expansions of successive orders in each element. The second one implies the determination of the most sensitive parameter and the third one imposes a minimum size of each element.
4.3 Sparse polynomial chaos
The efficiency of the sparse PC lies in the reduction of the number of PC coefficients that need to be evaluated. In this section, three different sparse PC expansions are presented, namely the sparse PC with lowrank index sets, the sparse PC with isotropic hyperbolic index sets and the sparse PC with anisotropic hyperbolic index sets. These methods are described in detail by Blatman [36]. Their major steps are recalled hereafter for the sake of clarity.
4.3.1 Lowrank index sets
In this approach, the polynomials ϕ _{ α } are retained or discarded from the PC expansion depending on their degrees p _{ α } and interaction orders j _{ α }, respectively defined as [25]$${p}_{\alpha}=\mathrm{}\left\right\alpha {}_{1}\text{\hspace{1em}}{j}_{\alpha}\mathrm{}={\displaystyle \sum _{i=1}^{r}}{1}_{\mathrm{}{\{\alpha}_{i}>0\}},\mathrm{}$$(22)where ${1}_{\{{\alpha}_{i}>0\}}$ is equal to 1 if α _{ i } > 0 and to 0 otherwise.
The number of retained PC coefficients is decreased by introducing a limited interaction order j ≤ p from which a reduced index set is defined [36]:$${\mathcal{A}}^{r,p,j}=\{\alpha \in {\mathbb{N}}^{\mathit{r}}:{p}_{\alpha}\le p,{j}_{\alpha}\le j\}\text{.}$$(23)
The lowrank PC expansion using the index set ${\mathcal{A}}^{\mathit{r},\mathit{p},\mathit{j}}$ is therefore defined as$${X}_{{\mathcal{A}}^{r,p,j}}(\xi )={\displaystyle \sum _{\alpha \in {\mathcal{A}}^{r,p,j}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}(\xi )\text{.}$$(24)
As explained in [36], the method is however inaccurate when j < min(r, p): in that case, the absence of the terms of interaction order greater than j prevents the PC expansion ${\mathit{X}}_{{\mathcal{A}}^{\mathit{r},\mathit{p},\mathit{j}}}$ from converging to the true value of the QoI. The sparse PC with hyperbolic index sets, presented in the next sections, overcomes this limitation.
4.3.2 Isotropic hyperbolic index sets
In this method, a reduced index set is introduced on the basis of a hyperbolic truncation using an mnorm with 0 < m < 1 (see [24])$${\mathcal{A}}_{m}^{r,p}=\{\alpha \in {\mathbb{N}}^{\mathit{r}}:\left\alpha \right{}_{m}\le p\}\text{,}$$(25)where$${\mathit{\alpha}}_{\mathit{m}}={\left({\displaystyle \sum _{i=1}^{r}}{\alpha}_{i}^{m}\right)}^{1/m}\text{.}$$(26)
Whatever the value of m chosen for the mnorm, the sequence of nested sets ${\mathcal{A}}_{\mathit{m}}^{\mathit{r},\mathit{p}}\text{\hspace{0.17em}}(\mathit{p}\in \mathbb{N})$ always converges to the set $\text{}{\mathbb{N}}^{\mathit{r}}$.
The isotropic hyperbolic PC expansion with the index set ${\mathcal{A}}_{\mathit{m}}^{\mathit{r},\mathit{p}}$ is then defined as$${X}_{{\mathcal{A}}_{m}^{r,p}}(\xi )={\displaystyle \sum _{\alpha \in {\mathcal{A}}_{m}^{r,p}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}(\xi )\text{.}$$(27)
4.3.3 Anisotropic hyperbolic index sets
In this third method, the strategy to truncate the PC expansions favors input random variables ξ _{ i } with large total sensitivity indices ${S}_{i}^{T}$. For this purpose, the truncation is based on the following anisotropic norm$$\left\right\alpha {}_{m,w}={\left({\displaystyle \sum _{i=1}^{r}}{w}_{i}{\alpha}_{i}{}^{m}\right)}^{1/m},\text{\hspace{1em}}{w}_{i}\ge 1\text{.}$$(28)
The corresponding anisotropic index set is then chosen as$${\mathcal{A}}_{m,w}^{r.p}=\{\alpha \in {\mathbb{N}}^{\mathit{r}}:\left\alpha \right{}_{m,w}\le p\}\text{,}$$(29)where w is a set of weights w _{ i } defined by$${w}_{i}=1+\frac{{\mathrm{max}}_{1\le j\le r}{S}_{j}^{T}{S}_{i}^{T}}{{\displaystyle \sum _{k=1}^{r}}{S}_{k}^{T}},\text{\hspace{1em}}i=\mathrm{1,\dots},r\text{.}$$(30)
In the above equation, ${S}_{i}^{T}$ is the PCbased total sensitivity index [37] of the QoI with respect to the input random variable ξ _{ i }, and is computed as$${S}_{i}^{T}=\frac{1}{{D}_{PC}}{\displaystyle \sum _{\alpha \in {I}_{i}^{+}}}{{\displaystyle \overline{x}}}_{\alpha}^{2}E[{\varphi}_{\alpha}^{2}(\xi )]\text{,}$$(31)where ${\mathcal{I}}_{i}^{+}$ denotes the set of indices having a nonzero ith component$${\mathcal{I}}_{i}^{+}=\left\{\alpha \in {\mathcal{A}}_{m,w}^{r,p}:{\alpha}_{i}\ne 0\right\}\text{,}$$(32)and D_{PC} the variance of the QoI$${D}_{PC}={\displaystyle \sum _{\alpha \in {\mathcal{A}}_{m,w}^{r,p}}}{{\displaystyle \overline{x}}}_{\alpha}^{2}E[{\varphi}_{\alpha}^{2}(\xi )]\text{.}$$(33)
The quantities ${S}_{i}^{T}$ approximately represent the amount of response variance that is explained by the variance of ξ _{ i }.
The anisotropic hyperbolic polynomial chaos expansions are finally defined with the index sets ${\mathcal{A}}_{m,w}^{r,p}$ as$${X}_{{\mathcal{A}}_{m,w}^{r,p}}(\xi )={\displaystyle \sum _{\alpha \in r{\mathcal{A}}_{m,w}^{r,p}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}(\xi )\text{.}$$(34)
4.3.4 Error estimates of the polynomial chaos approximations
The building of a sparse PC expansion is based on an incremental search of the significant terms, and therefore requires the use of error estimates to assess the accuracies of the consecutive PC approximations.
A theoretical error relevant in this context is the generalization error defined as$$\text{Err}=E[{(X(\xi ){{\displaystyle \widehat{\mathit{X}}}}_{\mathcal{A}}(\xi ))}^{2}]\text{,}$$(35)which is based on the difference between the deterministic evaluation X (ξ) of the QoI and its PC approximation ${\widehat{\mathit{X}}}_{\mathcal{A}}(\xi )$ computed from a finite non empty subset $\mathcal{A}\subset {\mathbb{N}}^{\mathit{r}}$, that is$${{\displaystyle \widehat{\mathit{X}}}}_{\mathcal{A}}(\xi )={\displaystyle \sum _{\alpha \in \mathcal{A}}}{{\displaystyle \overline{x}}}_{\alpha}{\varphi}_{\alpha}(\xi )\text{.}$$(36)
The generalization error is estimated in practice by the following empirical error:$${\text{Err}}_{\text{emp}}=\frac{1}{Q}{\displaystyle \sum _{q=1}^{Q}}\left[{(X({\xi}^{(q)}){\widehat{\mathit{X}}}_{\mathcal{A}}({\xi}^{(q)}))}^{2}\right]$$(37) in which the differences are computed specifically at the Q observations of an NED ${\xi}^{(q)}=\left({\xi}_{1}^{(q)},\mathrm{\dots},{\xi}_{r}^{(q)}\right)$. The latter will be used in the following to compute a coefficient of determination R ^{2} defined as$${R}^{2}=1\frac{{\text{Err}}_{\text{emp}}}{\widehat{{\displaystyle \mathit{V}}}[X]}\text{,}$$(38)where $\widehat{\mathit{V}}}[X]$ is the variance of X (ξ ^{ (q)}):$$\widehat{\mathit{V}}}[X]=\frac{1}{Q1}{\displaystyle \sum _{q=1}^{Q}}{(X({\xi}^{(q)})\overline{\mathit{X}})}^{2}\text{\hspace{1em}}\text{with}\text{\hspace{1em}}{\displaystyle \overline{X}}=\frac{1}{Q}{\displaystyle \sum _{q=1}^{Q}}X({\xi}^{(q)})\text{.$$
An overfitting phenomenon is known to occur when using the empirical error, which, as a consequence, underestimates the generalization error. The leaveoneout error [36], which is based on a sum of squared predicted residuals Δ^{(i)} defined hereafter, may be useful to avoid this drawback. A predicted residual expresses the difference between the deterministic evaluation X(ξ ^{(i)}) of the QoI at the ith observation of the NED ξ ^{(q)}, and its prediction $\widehat{\mathit{X}}{\mathrm{}}_{\mathcal{A}}^{(\mathrm{i})}({\xi}^{(\mathrm{}i)\mathrm{}})\mathrm{}$ obtained with a PC expansion $\widehat{\mathit{X}}{\mathrm{}}_{\mathcal{A}}^{(\mathrm{i})}\mathrm{}$ built from a reduced NED (ξ ^{(1)},…, ξ ^{(Q)}) ξ ^{(i)} (that is the original NED from which the observation ξ ^{ (i)} has been discarded) [25]:$${\mathrm{\Delta}}^{\mathrm{}(i\mathrm{})}=\mathrm{}X({\xi}^{\mathrm{}(i)\mathrm{}})\mathrm{}{{\displaystyle \widehat{\mathit{X}}}}_{\mathcal{A}}^{\mathrm{}(i)\mathrm{}}(\mathrm{}{\xi}^{\mathrm{}(i\mathrm{})})\mathrm{}\text{.}$$(39)
The leaveoneout error is then defined as$${\text{Err}}_{\text{LOO}}=\frac{1}{Q}{\displaystyle \sum _{i=1}^{Q}}{({\mathrm{\Delta}}^{(i)})}^{2}\text{.}$$(40)
In practice, the predicted residual ${\mathrm{\Delta}}^{(\mathrm{}i)\mathrm{}}$ may be computed as [36]$${\mathrm{\Delta}}^{(i)}=\frac{X({\xi}^{(i)}){\widehat{\mathit{X}}}_{\mathcal{A}}({\xi}^{(i)})}{1{h}_{i}}\text{,}$$(41)where h _{ i } is the ith diagonal term of the matrix $\varphi ({\xi}^{(q)}){({\varphi}^{T}({\xi}^{(q)})\varphi ({\xi}^{(q)}))}^{1}{\varphi}^{T}({\xi}^{(q)})$. The leaveoneout error is in that case given by$${\text{Err}}_{\text{LOO}}=\frac{1}{Q}{\displaystyle \sum _{i=1}^{Q}}{\left(\frac{X({\xi}^{(i)}){\widehat{\mathit{X}}}_{\mathcal{A}}({\xi}^{(i)})}{1{h}_{i}}\right)}^{2}\text{.}$$(42)
A determination coefficient S ^{2} equivalent to that of the empirical error, R ^{2}, may be defined for the leaveonout error:$${S}^{2}=1\frac{{\text{Err}}_{\text{LOO}}}{\widehat{\mathit{V}}[X]}\text{.}$$(43)
The two coefficients R ^{2} and S ^{2} defined above will be used in an algorithm whose aim is to build an optimal sparse PC expansion involving the most significant terms from an adapted NED of reduced size. This algorithm is described in the next section.
4.3.5 Sparse PC expansion building algorithm
As explained previously, the efficiency of the method may be increased by retaining only the most significant PC polynomials [25] among those corresponding to the index sets ${\mathcal{A}}^{r,p,j}$, ${\mathcal{A}}_{m}^{r,p}$ and ${\mathcal{A}}_{m,w}^{r,p}$. In the following, the final index sets of the kept terms are respectively denoted as ${\mathcal{A}}^{p}$, ${\mathcal{A}}_{m}^{p}$ and ${\mathcal{A}}_{m,w}^{p}$.
The search for those most significant coefficients is performed through an iterative procedure which is summarized below in 5 basic steps.
Step 1

Select an NED (ξ ^{ (q)}), e.g. a random design based on LHS [33], of arbitrary size Q _{ k } = 4N _{ p }, where N _{ p } is determined by equation (15) with r uncertain parameters and p = 1. The DM evaluations at the NED points are gathered in the vector X (ξ ^{ (q)}).

Set arbitrarily the values of the parameters corresponding to the chosen sparse PC method: the maximal PC order p _{max} and interaction order j _{max} with lowrank index sets (or the coefficient m used for the mnorm of truncation for isotropic and anisotropic hyperbolic index sets), as well as the target accuracy ${S}_{\text{target}}^{2}$ and two thresholds ϵ _{1} and ϵ _{2}.
Step 2
Initialize the algorithm: the PC order is set to p = 0, and the truncation index set to the null element of ${\mathbb{N}}^{\mathit{r}}$, {0} (the vector of weights w _{ i } is set to w = {1,…, 1} for anisotropic index steps). The corresponding initial values of the determination coefficients are denoted as ${R}_{0}^{2}$ and ${S}_{0}^{2}$.
Step 3: Training step − enrichment of the PC basis
Increment the order value: p → p + 1 ∈ [1,…, p _{max}].
$\Rightarrow \mathbf{F}\mathbf{o}\mathbf{r}\mathbf{w}\mathbf{a}\mathbf{r}\mathbf{d}\text{\hspace{0.17em}}\mathbf{s}\mathbf{t}\mathbf{e}\mathbf{p}$ (addition step):

With lowrank index sets: set the interaction order value to $j\mathrm{}=\mathrm{min}(\mathrm{}p\mathrm{},\text{\hspace{0.17em}}{j}_{\text{max}})\mathrm{}$; for each term from the candidate set ${{\mathcal{C}}^{\mathit{j},\mathit{p}}=\{\mathit{\alpha}\in \mathbb{N}}^{\mathit{r}}:{\mathit{p}}_{\mathit{\alpha}}=\mathit{p},{\mathit{j}}_{\mathrm{\alpha}}=\mathit{j}\}\text{\hspace{0.17em}}$ , add it to the set ${\mathcal{A}}^{\mathit{p}1}$, compute the PC coefficients by regression (Eq. (17)) using the Q _{ k } points and the corresponding determination coefficient R ^{2}.
With isotropic (resp. anisotropic) index sets: for each term from the candidate set ${\mathcal{C}}_{m}^{p}=\{\alpha \in {\mathbb{N}}^{r}:p1\le \left\alpha \right{}_{m}\le p\}$ (resp. ${\mathcal{C}}_{m,w}^{p}=\{\alpha \in {\mathbb{N}}^{r}:p1\le \left\alpha \right{}_{m,w}\le p\}$), add it to the set ${\mathcal{A}}_{m}^{p1}$ (resp. ${\mathcal{A}}_{m,w}^{p1}$) and compute, as above, the PC coefficients and the determination coefficient R ^{2}.

Retain each candidate whose addition leads to a significant increase in the value of the coefficient R ^{2}, that is$$\mathrm{\Delta}{R}^{2}={R}^{2}{R}_{0}^{2}\ge {\u03f5}_{1}\text{.}$$(44)
If inequality (44) is not respected the candidate is discarded.
With lowrank index sets, let ${\mathcal{A}}^{p,+}$ be the final truncation set at this stage (respectively ${\mathcal{A}}_{m}^{p+}$ and ${\mathcal{A}}_{m,w}^{p+}$ for isotropic and anisotropic index sets).
⇒ Backward step (elimination step):

Remove in turn each term in ${\mathcal{A}}^{p,+}$ (${\mathcal{A}}_{m}^{p+}$ and ${\mathcal{A}}_{m,w}^{p+}$ respectively with isotropic and anisotropic index sets) of order strictly lower than p, and compute again the PC expansion coefficients and the associated coefficient R ^{2} in each case.

Discard from ${\mathcal{A}}^{p,+}$ (${\mathcal{A}}_{m}^{p+}$ and ${\mathcal{A}}_{m,w}^{p+}$ respectively with isotropic and anisotropic index sets) the terms that lead to an insignificant decrease in R ^{2}, i.e. if the suppression of the candidate leads to$$\mathrm{\Delta}{R}^{2}={R}_{0}^{2}{R}^{2}<{\u03f5}_{2}\text{.}$$(45)
If inequality (45) is not respected the candidate is kept.
Let ${\mathcal{A}}^{p}$ be the final truncation set for lowrank index sets (respectively ${\mathcal{A}}_{m}^{p}$ and ${\mathcal{A}}_{m,w}^{p}$ for isotropic and anisotropic index sets). With anisotropic index sets, the total sensitivity indices ${S}_{i}^{T}$ of the current PC approximation are computed and the weights w _{ i } are updated (Eq. (30)).
Step 4: Verification of the conditioning of the regression information matrix
1. If the conditioning is satisfying, i.e. the size Q _{ k } of the NED (ξ ^{ (q)}) is larger than 2.card(${\mathcal{A}}^{p}$) with lowrank index sets (respectively 2.card(${\mathcal{A}}_{m}^{p}$) and 2.card(${\mathcal{A}}_{m,w}^{p}$) for isotropic and anisotropic index sets) go to step 5.
2. If the conditioning is poor, i.e. the size Q _{ k } of the NED (ξ ^{ (q)}) is smaller than 2.card(${\mathcal{A}}^{p}$) with lowrank index sets (respectively 2.card(${\mathcal{A}}_{m}^{p}$) and 2.card(${\mathcal{A}}_{m,w}^{p}$) for isotropic and anisotropic index sets), an enrichment of the NED is done using nested Latin Hypercube designs [25,38] to reach a size Q _{ k+1}. In this case, the truncation set is reset to {0} and the enrichment procedure is restarted from step 2.
Step 5: Test step
Stop if either the leaveoneout error ${S}_{0}^{2}$ is larger than the target value ${S}_{\mathrm{t}\mathrm{arg}\mathrm{e}\mathrm{t}}^{2}$ or if the order of the PC expansion is equal to p _{max}. Otherwise, go back to step 3.
The detailed algorithms are presented in Appendix A.
5 Application and results
In this section, the stability of the clutch system with uncertain parameters is studied using the lumped model presented in Section 2.1. The analysis is carried out using MCDM as the reference solution and MC computations coupled with polynomial chaos methods (i.e. gPC, MEgPC and sparse PC). The comparisons between the gPC expansions and the reference results is performed in terms of statistics (mean and variance) of the propensity of stability, which is defined hereafter.
5.1 Propensity of stability
The propensity of stability V of the clutch system is estimated for a given random space of uncertain parameters with the MC method from the following equation:$$V=\frac{{N}_{\text{stable}}}{N}\text{,}$$(46)where N is the total number of samples and N _{stable} the number of stable samples.
To ensure a confidence level of 99% with an accuracy margin of 1.07% for the stability propensity, a total number of samples N = 10 000 is required for a given space of the stochastic parameters [39]. It should be noted that this sample number N does not depend on the number of uncertain parameters, but only on the desired confidence level and accuracy margin.
If the computational cost to build the polynomial chaos is greater than 10 000 Calculations with the DM (CDM), there is no advantage in using the developments in polynomial chaos. In this case, MCDM remains the optimal method. The limiting computational cost for building chaos is therefore 10 000 CDM.
5.2 Uncertain parameters
In the analytical model of the clutch system with 6 DOF presented in Section 2.1, eight parameters may be chosen uncertain: μ, k _{ p }, k _{ f }, γ _{1}, γ _{2}, r _{1}, r _{2} and l. The uncertain parameters are considered independent, random and uniform within the intervals [0.95V _{ n }, 1.05V _{ n }] where V _{ n } refers to the nominal values of the parameters (given at the end of Sect. 2.1); μ is supposed independent, random and uniform in the interval [0, 0.5]. Computations of the clutch system eigenvalues have been performed in eight cases corresponding to varying numbers of uncertain parameters r = 1 to 8. The lists of uncertain parameters for the different cases are given in Table 2, the other parameters being set to their nominal values.
Mean and variance of the real part of λ _{2} and propensity of stability obtained using the MCDM.
5.3 Reference study: Monte Carlo on the deterministic model
In this section, the MCDM is applied to study the stability of the clutch system. To this end, 10 000 sets of eigenvalues λ _{ i } (i = 1,…, 12), solutions of the characteristic equation (9), are directly calculated from 10 000 independent random samples.
After computation, it appears that the clutch system has 8 nonzero eigenvalues (modes 1, 2, 3 and 4 with complex conjugates). Analyzing the eigenvectors corresponding to the 4 quasizero eigenvalues we could see that these modes are mostly constituted of wobbling and nodaldiameter bending motions. Modes 3 and 4 are seen always decoupled and stable, contrary to modes 1 and 2 which are affected by a coalescence phenomenon and therefore rule the stability of the clutch system.
Figure 2 displays the real and imaginary parts of modes 1 (gray) and 2 (black) computed with 10 000 values of the friction coefficient μ, the other parameters being set to their nominal values (that is r = 1). When μ increases from 0 to 0.105, the imaginary parts (i.e. 2π× frequencies) of both eigenvalues are separated but tend to come closer, while the real parts are equal and negative, meaning that the system is stable. From μ = 0.105, the imaginary parts of the two modes overlap, characterizing the coalescence phenomenon, while the real parts separate. Furthermore, from the Hopf bifurcation point (μ = 0.11), the real part of mode 2 becomes positive and the system is unstable.
It is wellknown from the literature that for friction systems, the coalescence patterns may be different according to the values of the damping coefficients [5,40]. Two main effects may occur: a shift of the real parts of the eigenvalues called the shifting (or lowering) effect, and a smoothing effect in the vicinity of the coalescence point. The first one stabilizes the system whereas the second one can lead to an unintuitive destabilization of the system.
Trinh et al. have studied these effects in the case of the clutch system model used in this paper [20], for varying values of the damping coefficients. Their appearances proved to be strongly dependent on the ratio between the modal damping coefficients of the two coalescent modes.
The damping coefficients chosen in the present paper induce very close values of those two modal damping coefficients. Therefore, only the lowering effect appears. Small variations of the damping values due to parameter uncertainties may have very little influence on the coalescence patterns. However, the latter remain similar to those shown in Figure 2, and cannot affect the PC expansion efficiency. The influence of the lowering and smoothing effects on the efficiency of the different sparse PC expansions could be studied in a future paper. Modeling the eigenvalues from sparse PC expansions is expected to be easier if the smoothing effect is predominant, as the evolutions of their real parts are smoother in that case.
Eight MCDM studies have been performed to compute the QoI (i.e. the system eigenvalues) for the different cases of uncertain parameters (r = 1 to 8). The resulting statistics (mean, variance) of the real part of λ _{2} along with the corresponding propensity of stability are presented in Table 2.
Fig. 2 (a) Real parts and (b) imaginary parts of eigenvalues λ _{1} and λ _{2} for modes 1 and 2. 
5.4 Sparse PC method application for the stability analysis of a clutch system
In this section, the stability analysis is performed by coupling the MC method with the sparse PC expansions detailed in Section 4. Results obtained previously with classic gPC and MEgPC on the same model [20] are also presented for comparison purposes.
More specifically, the eigenvalues λ _{ i } (i = 1,…, 8) are first expanded into the following form$${\lambda}_{i}\mathrm{}(\xi \mathrm{})\approx {\displaystyle \sum _{\alpha \in \mathcal{A}}}{{\displaystyle \overline{\lambda}}}_{i,\mathrm{}\alpha}{\varphi}_{\alpha}\mathrm{}(\xi )\mathrm{}$$(47)where $\mathcal{A}$ denotes identically ${\mathcal{A}}^{\mathit{r},p}$ (for gPC and MEgPC), ${\mathcal{A}}^{p}$ (for sparse PC with low rank index sets), ${\mathcal{A}}_{m}^{p}$ or ${\mathcal{A}}_{m,w}^{p}$ (for isotropic and anisotropic hyperbolic index sets, respectively).
In a second step, the statistics (such as the mean and the variance of the real parts of the eigenvalues) are obtained within the MC framework by computing the eigenvalues using equation (47) for 10 000 random samples ξ.
Table 3 compares the results (in terms of stability propensity) obtained with MCDM, gPC, MEgPC, and the three sparse PC methods that are respectively referred to as SgPCLR (sparse PC with lowrank index sets), SgPCIH (sparse PC with isotropic hyperbolic index sets) and SgPCAH (sparse PC with anisotropic hyperbolic index sets).
The results obtained with gPC and MEgPC have been computed by Trinh et al. [20] using the regression method with Q = (p + 1)^{ r }. In that study, the PC expansions were built within an iterative procedure and the choice of the order p was based on a criterion which considered the decay rate of the relative error between two developments of successive orders.
Regarding gPC, the above strategy leads to a relatively low order p which starts from p = 6 for r = 1 and decreases to p = 2 for higher numbers of uncertain parameters, as shown in Table 3. The number of polynomials in the expansions become nevertheless excessive from r = 5, leading to a prohibitive number of CDM to evaluate the corresponding PC coefficients. The relative errors of the stability propensity between gPC and MCDM become simultaneously higher than 20%, the gPC being therefore inefficient in these cases.
The MEgPC method gives satisfying results for r ≤ 7 regarding the propensity of stability. However, with more than 5000 CDM required, the computational cost becomes expensive for r = 5 →7. For r = 8, the MEgPC is not worth using anymore as the number of CDM to compute the PC coefficients exceeds the limit of 10 000.
In the present study, sparse PC expansions have been computed with a target accuracy ${S}_{\text{target}}^{2}=0.999$, a maximal PC order p _{max} = 15 and two identical threshold values ${\mathit{\epsilon}}_{1}\mathrm{}={\mathit{\epsilon}}_{2}=0.001(1{S}_{\text{target}}^{2}\mathrm{})$. For the SgPCLR (with lowrank index sets), the following maximal interaction orders have been used: j _{max} = 2 for r = 3, j _{max} = 3 for r = 4, 5, j _{max} = 4 for r = 6, 7 and j _{max} = 5 for r = 8. When computing sparse PC expansions with hyperbolic index sets, two different values of the parameter m defining the mnorm have been used: in the isotropic case (SgPCIH), the value of m was set to 0.6, a higher value leading to an excessive number of CDM; in the anisotropic case (SgPCAH), a value of m = 0.8 enabled to limit the errors of the stability propensity while keeping a reasonable number of CDM. Table 3 shows the results obtained with the three kinds of sparse PC expansions for r = 3 →8, as those methods are not relevant for a lower number of uncertain parameters. Moreover, in order to evaluate the influence of the choice of the initial NED on the results, expansions from 5 distinct NED have been computed for each sparse PC method. Table 3 only displays in square brackets, in each case, the minimal and maximal values of the PC orders, the CDM numbers and the relative errors.
In this study, the sparse methods are seen highly efficient compared to the gPC and MEgPC, especially for r ≥ 5, as they ensure a high accuracy while using a smaller number of CDM. All the relative errors of the stability propensity between sparse PC and MCDM are lower than 5% for r = 3 →5. For r = 6 →8 the relative errors using SgPCLR and SgPCAH are still lower than 5%, and lower than 7% for SgPCIH. For r = 8, the efficiency of the sparse methods is particularly visible as the number of CDM exceeds the limit of 10 000 with MEgPC, while the relative error becomes excessive with gPC. With the SgPCLR, for most initial NED the number of CDM is lower than 5000, meaning a reduction by about 50% of the computational cost, and the accuracy reaches a high level with relative errors lower than 1%. The number of CDM using SgPCIH and SgPCAH are lower than 3000, leading to a computational cost reduction of more than 70%; the relative errors are fully satisfying but slightly higher than with SgPCLR (<5% with SgPCAH and <7% with SgPCAH). Although the global accuracy of the results with sparse PC expansions is high whatever the value of r, the choice of NED is seen to have a significant influence on the resulting number of CDM and stability propensity values.
Comparing the different sparse PC methods to one another, it appears that the results are the most accurate with the SgPCLR method whatever the number of uncertain parameters. However, the computational cost increases substantially with this method compared to SgPCIH and SgPCAH. Between the two sparse PC methods with hyperbolic index sets, i.e. SgPCIH and SgPCAH, the accuracy is seen significantly improved with the SgPCAH for similar computational costs. This behavior is favored by the use of weights for the uncertain parameters. The total sensitivity indices ${S}_{i}^{T}$ and the weights w _{ i } used for SgPCAH are shown in Table 4. As previously, the results shown (in square brackets) are restricted to the minimal and maximal values obtained for the 5 NED. The table shows that the friction coefficient μ is the most sensitive parameter; indeed, its sensitivity index is at least 400 times greater than the sensitivity indices of the other parameters. Thus, its weighting is the smallest (1 instead of about 1.95). The higher the total sensitivity indices of uncertain parameters are, the smaller are their weights, which increases the influence of those important parameters on the PC coefficients, and therefore the accuracy of the results, particularly with a higher number of uncertain parameters.
To illustrate the high accuracy of the results obtained with sparse PC methods, Figure 3 shows the real parts of the eigenvalues corresponding to modes 1 and 2 using MCDM (black) and SgPCAH (gray) for r = 3 to 8. The results displayed correspond in each case to the initial NED that leads to the smallest relative error regarding the stability propensity. The details of the corresponding expansions (order p, number of CDM and relative error) are specified in Table 5. Each SgPCAH expansion is computed for r uncertain parameters, but in Figure 3 the QoI are plotted with respect to the sole friction coefficient μ for given values of the r − 1 other parameters. The plots show a good agreement between the SgPCAH and the reference solutions. In particular, we can see that the Hopf bifurcation point is well determined by the sparse PC approximation, leading to an accurate estimation of the propensity of stability.
Comparison of the results between gPC, MEgPC, sparse PC and MCDM.
Total sensitivity indices and weights of uncertain parameters with SgPCAH.
Parameters of the SgPCAH expansions.
Fig. 3 Real parts of the eigenvalues of modes 1 and 2 using MCDM (black) and SgPCAH (gray) as functions of the friction coefficient μ. 
6 Conclusion
In this paper, several sparse PC expansions have been used to analyze the stability of a 6 DOF clutch system. The sparse PC methods are compared to the classic MC method performed directly on the Deterministic clutch Model (MCDM) and to usual gPC and MEgPC.
To overcome the limitations of the classic MC approach which is often too costly for large industrial systems, gPC and MEgPC have been used in the past; however, the computational costs with these methods become also prohibitive when the number of uncertainties is large (r ≥ 5). In this study, three sparse PC expansions are proposed to reduce those costs, notably by limiting the number of polynomials which express an interaction between several uncertain parameters. In each case, the sparse PC expansion is built using an iterative procedure based on error computations which retains only the most significant coefficients among the initial reduced index set of the polynomials, and adjusts automatically the size of the NED used to evaluate them. The resulting final set of PC coefficients, as well as the number of simulations of the DM required to evaluate them, are consequently extremely reduced in comparison with full PC approximations.
The results show that the use of sparse PC methods to analyze the stability of a clutch system with uncertain parameters is efficient for a number of uncertain parameters ranging from 3 to 8. In these cases, it allows a substantial reduction of the computational cost while ensuring a high accuracy of the results in comparison with the MCDM, gPC and MEgPC methods. The numerous tests carried out in this study however show that the choice of the initial NED has a noticeable influence on the results, such as the number of required simulations of the DM, which can be more than doubled between the worst and the most favorable configurations in the case of the sparse PC method with lowrank index sets. The benefits associated with the use of a sparse PC expansion yet remain interesting in all the cases. Given those promising results, these methods could also be used advantageously to study dynamical systems with numerous DOF and uncertain parameters.
Appendix A Algorithms applied with sparse PC
The different algorithms used to develop the sparse PC expansions are detailed respectively in Figure A1 for lowrank index sets, in Figure A2 for isotropic hyperbolic index sets and in Figure A3 for anisotropic hyperbolic index sets.
Fig. A1 Algorithm applied to build a sparse polynomial chaos expansion with lowrank index sets. 
Fig. A2 Algorithm applied to build a sparse polynomial chaos expansion with isotropic hyperbolic index sets. 
Fig. A3 Algorithm applied to build a sparse polynomial chaos expansion with anisotropic hyperbolic index sets. 
References
 B. Hervé, J.J. Sinou, H. Mahé, L. Jézéquel, Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear selfexcited system subject to gyroscopic and circulatory actions, J. Sound Vib. 323 (2009) 944–973 [Google Scholar]
 F. Van De Velde, P. De Baets, A new approach of stickslip based on quasiharmonic tangential oscillations, Wear 216 (1998) 15–26 [Google Scholar]
 F. Van De Velde, P. De Baets, The relation between friction force and relative speed during the slipphase of a stickslip cycle, Wear 219 (1998) 220–226 [Google Scholar]
 P. Wickramarachi, R. Singh, G. Bailey, Analysis of frictioninduced vibration leading to “eek” noise in a dry friction clutch, Noise Control Eng. J. 53 (2005) 138–144 [CrossRef] [Google Scholar]
 G. Fritz, J.J. Sinou, J.M. Duffal, L. Jézéquel, Investigation of the relationship between damping and modecoupling patterns in case of brake squeal, J. Sound Vib. 307 (2007) 591–609 [Google Scholar]
 B. Hervé, J.J. Sinou, H. Mahé, L. Jézéquel, Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects, Eur. J. Mech. A Solids 27 (2008) 141–160 [Google Scholar]
 D. Centea, H. Rahnejat, M.T. Menday, Nonlinear multibody dynamic analysis for the study of clutch torsional vibrations (judder), Appl. Math. Model. 25 (2001) 177–192 [Google Scholar]
 G. Chevallier, F. Macewko, F. RobbeValloire, Dynamic friction evolution during transient sliding, Tribol. Ser. 43 (2003) 537–543 [CrossRef] [Google Scholar]
 G. Chevallier, D. Le Nizerhy, Friction induced vibrations in a clutch system. consequences on the apparent friction torque, in: C.A. Motasoares, J.A.C. Martins, H.C. Rodrigues, J.A.C. Ambrósio, C.A.B. Pina, C.M. Motasoares, E.B.R. Pereira, J. Folgado (Eds.), III European Conference on Computational Mechanics: Solids, Structures and Coupled Problems in Engineering: Book of Abstracts, Dordrecht, Springer, Netherlands, 2006, 309–309 [Google Scholar]
 R.A. Ibrahim, Frictioninduced vibration, chatter, squeal, and chaos − Part I: mechanics of contact and friction, Appl. Mech. Rev. 47 (1994) 209–226 [Google Scholar]
 L. Nechak, S. Berger, E. Aubry, A polynomial chaos approach to the robust analysis of the dynamic behaviour of friction systems, Eur. J. Mech. A Solids 30 (2011) 594–607 [Google Scholar]
 L. Nechak, S. Berger, E. Aubry, Nonintrusive generalized polynomial chaos for the robust stability analysis of uncertain nonlinear dynamic friction systems, J. Sound Vib. 332 (2013) 1204–1215 [Google Scholar]
 E. Sarrouy, O. Dessombz, J.J. Sinou, Piecewise polynomial chaos expansion with an application to brake squeal of a linear brake system, J. Sound Vib. 332 (2013) 577–594 [Google Scholar]
 E. Sarrouy, O. Dessombz, J.J. Sinou, Stochastic analysis of the eigenvalue problem for mechanical systems using polynomial chaos expansion − application to a finite element rotor, J. Vib. Acoust. 134 (2012) 051009 [Google Scholar]
 L. Nechak, S. Besset, J.J. Sinou, Robustness of stochastic expansions for the stability of uncertain nonlinear dynamical systems − application to brake squeal, Mech. Syst. Signal Process. 111 (2018) 194–209 [Google Scholar]
 L. Nechak, J.J. Sinou, Hybrid surrogate model for the prediction of uncertain frictioninduced instabilities, J. Sound Vib. 396 (2017) 122–143 [Google Scholar]
 E. Denimal, L. Nechak, J.J. Sinou, S. Nacivet, A novel hybrid surrogate model and its application on a mechanical system subjected to frictioninduced vibration, J. Sound Vib. 434 (2018) 456–474 [Google Scholar]
 E. Jacquelin, O. Dessombz, J.J. Sinou, S. Adhikari, M.I. Friswell, Polynomial chaosbased extended padé expansion in structural dynamics, Int. J. Numer. Methods Eng. 111 (2017) 1170–1191 [Google Scholar]
 L. Nechak, S. Berger, E. Aubry, Prediction of random self frictioninduced vibrations in uncertain dry friction systems using a multielement generalized polynomial chaos approach, J. Vib. Acoust. 134 (2012) 041015 [Google Scholar]
 M.H. Trinh, S. Berger, E. Aubry, Stability analysis of a clutch system with multielement generalized polynomial chaos, Mech. Ind. 17 (2016) 205 [CrossRef] [Google Scholar]
 D.C. Montgomery, Design and analysis of experiments, in: Student solutions manual, John Wiley & Sons, NY, 2008 [Google Scholar]
 J. An, A. Owen, Quasiregression, J. Complexity 17 (2001) 588–607 [CrossRef] [Google Scholar]
 R.A. Todor, C. Schwab, Convergence rates for sparse chaos approximations of elliptic problems with stochastic coefficients, IMA J. Numer. Anal. 27 (2007) 232–261 [CrossRef] [MathSciNet] [Google Scholar]
 G. Blatman, B. Sudret, Anisotropic parcimonious polynomial chaos expansions based on the sparsityofeffects principle, in: Proc ICOSSAR'09, International Conference in Structural Safety and Relability, 2009 [Google Scholar]
 G. Blatman, B. Sudret, An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probabilist. Eng. Mech. 25 (2010) 183–197 [Google Scholar]
 G. Blatman, B. Sudret, Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach, Comptes Rendus Mécanique 336 (2008) 518–523 [CrossRef] [Google Scholar]
 J.J. Sinou, Transient nonlinear dynamic analysis of automotive disc brake squeal − On the need to consider both stability and nonlinear analysis, Mech. Res. Commun. 37 (2010) 96–105 [Google Scholar]
 A. Shapiro, Monte Carlo sampling methods, in: Stochastic Programming, Vol. 10, Handbooks in Operations Research and Management Science, Elsevier Science, NY, 2003, 353–425 [Google Scholar]
 D. Xiu, G.E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2002) 619–644 [Google Scholar]
 R.H. Cameron, W.T. Martin, The orthogonal development of nonlinear functionals in series of fourierhermite functionals, Ann. Math. 48 (1947) 385–392 [Google Scholar]
 R. Askey, J.A. Wilson, Some basic hypergeometric orthogonal polynomials that generalize Jacobi polynomials, Mem. Am. Math. Soc. 319 (1985) [Google Scholar]
 D. Xiu, G.E. Karniadakis, Modeling uncertainty in steady state diffusion problems via generalized polynomial chaos, Comput. Methods Appl. Mech. Eng. 191 (2002) 4927–4948 [CrossRef] [MathSciNet] [Google Scholar]
 M.D. McKay, R.J. Beckman, W.J. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21 (1979) 239–245 [Google Scholar]
 M. Berveiller, B. Sudret, M. Lemaire, Stochastic finite element: a non intrusive approach by regression, Eur. J. Comput. Mech. (Revue Européenne de Mécanique Numérique) 15 (2006) 81–92 [CrossRef] [Google Scholar]
 X. Wan, G.E. Karniadakis, An adaptive multielement generalized polynomial chaos method for stochastic differential equations, J. Comput. Phys. 209 (2005) 617–642 [Google Scholar]
 G. Blatman, Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis, PhD thesis, Université BLAISE PASCAL − Clermont II, 2009 [Google Scholar]
 B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliab. Eng. Syst. Saf. 93 (2008) 964–979 [CrossRef] [Google Scholar]
 G.G. Wang, Adaptive response surface method using inherited Latin Hypercube design points, J. Mech. Des. 125 (2003) 210–220 [CrossRef] [Google Scholar]
 G. Baillargeon, Méthodes statistiques de l'ingénieur, Les Editions SMG, 1990 [Google Scholar]
 N. Hoffmann, L. Gaul, Effects of damping on modecoupling instability in friction induced oscillations, J. Appl. Math. Mech. 83 (2003) 524–534 [Google Scholar]
Cite this article as: D.T. Kieu, B. Bergeot, M.L. Gobert, S. Berger, Stability analysis of a clutch system with uncertain parameters using sparse polynomial chaos expansions, Mechanics & Industry 20, 104 (2019)
All Tables
Correspondence between the families of orthogonal polynomials and the distributions of the random variables.
Mean and variance of the real part of λ _{2} and propensity of stability obtained using the MCDM.
All Figures
Fig. 1 Analytical model of the clutch system. On the left two rotations θ _{ x }, θ _{ y } (blue arrows) around the x and y axes, on the right four translation movements Z _{ A }, Z _{ B }, Z _{ C }, Z _{ D } (green arrows) of points A, B, C, D in the z direction. 

In the text 
Fig. 2 (a) Real parts and (b) imaginary parts of eigenvalues λ _{1} and λ _{2} for modes 1 and 2. 

In the text 
Fig. 3 Real parts of the eigenvalues of modes 1 and 2 using MCDM (black) and SgPCAH (gray) as functions of the friction coefficient μ. 

In the text 
Fig. A1 Algorithm applied to build a sparse polynomial chaos expansion with lowrank index sets. 

In the text 
Fig. A2 Algorithm applied to build a sparse polynomial chaos expansion with isotropic hyperbolic index sets. 

In the text 
Fig. A3 Algorithm applied to build a sparse polynomial chaos expansion with anisotropic hyperbolic index sets. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.