© D.T. Kieu et al., published by EDP Sciences 2019

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 stick-slip 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 sprag-slip or to the mode-coupling 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 stick-slip [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 multi-element generalized polynomial chaos (ME-gPC), 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 non-intrusive methods. The intrusive approach requires to modify the DM to determine the PC coefficients. In the non-intrusive 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 non-intrusive 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 non-intrusive 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 non-intrusive gPC with respectively the inclusion function based on Chebyshev polynomials and the kriging formalism are considered to predict friction-induced 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 ME-gPC may be used as an alternative. For instance, the ME-gPC appeared efficient to predict the friction-induced 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 ME-gPC 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 ME-gPC expansions, may be considered. They are inspired by the so-called sparsity-of-effects principle [21], which states that most models are principally governed by main effects and low-order 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 low-rank 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 low-rank 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 ME-gPC 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, ME-gPC 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 mode-coupling 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 γ ₁; a second factor γ ₂ is similarly used at locations Cʹ and D ʹ, as detailed in equation (1) below [4]$$k_A={\gamma }_1{k}_p/4,\text{}{k}_B=k_p/(4{\gamma }_1)\text{,}$$(1a) $$k_C={\gamma }_2{k}_p/4,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 ₁ and external radius r ₂. The points Aʹ, B ʹ, C ʹ and D ʹ are supposed located at the average radius ${\displaystyle \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 rigid-body 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 nodal-diameter 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),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),N_{D^{\prime }}=k_D(Z_D-{\displaystyle \overline{r}}{\theta }_y)\text{,}$$(2b)and$$F_{A^{\prime }}=\mu N_{A^{\prime }},F_{B^{\prime }}=\mu N_{B^{\prime }}\text{,}$$(3a) $$F_{C^{\prime }}=\mu N_{C^{\prime }},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}}{\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+4k_f\right)\hfill & \hfill \mu l{\displaystyle \bar{r}}\left(k_C+k_D\right)\hfill & \hfill {\displaystyle \bar{r}}\left(k_A+2k_f\right)\hfill & \hfill -{\displaystyle \bar{r}}\left(k_B+2k_f\right)\hfill & \hfill \mu l{k}_C\hfill & \hfill -\mu l{k}_D\hfill \\ \hfill -\mu l{\displaystyle \bar{r}}\left(k_A+k_B\right)\hfill & \hfill {{\displaystyle \bar{r}}}^2\left(k_C+k_D+4k_f\right)\hfill & \hfill -\mu l{k}_A\hfill & \hfill \mu l{k}_B\hfill & \hfill {\displaystyle \bar{r}}\left(k_C+2k_f\right)\hfill & \hfill -{\displaystyle \bar{r}}\left(k_D+2k_f\right)\hfill \\ \hfill {\displaystyle \bar{r}}(k_A+2k_f)\hfill & \hfill 0\hfill & \hfill k_A+2k_f\hfill & \hfill 0\hfill & \hfill -k_f\hfill & \hfill -k_f\hfill \\ \hfill -{\displaystyle \bar{r}}(k_B+2k_f)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_B+2k_f\hfill & \hfill -k_f\hfill & \hfill -k_f\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \bar{r}}(k_C+2k_f)\hfill & \hfill -k_f\hfill & \hfill -k_f\hfill & \hfill k_C+2k_f\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\displaystyle \bar{r}}(k_D+2k_f)\hfill & \hfill -k_f\hfill & \hfill -k_f\hfill & \hfill 0\hfill & \hfill k_D+2k_f\hfill \end{array}\right]$$(7)

and $$\text{C}=\left[\begin{array}{cccccc}\hfill {{\displaystyle \bar{r}}}^2\left(c_A+c_B+4c_f\right)\hfill & \hfill \mu l{\displaystyle \bar{r}}\left(c_C+c_D\right)\hfill & \hfill {\displaystyle \bar{r}}\left(c_A+2c_f\right)\hfill & \hfill -{\displaystyle \bar{r}}\left(c_B+2c_f\right)\hfill & \hfill \mu l{c}_C\hfill & \hfill -\mu l{c}_D\hfill \\ \hfill -\mu l{\displaystyle \bar{r}}\left(c_A+c_B\right)\hfill & \hfill {{\displaystyle \bar{r}}}^2\left(c_C+c_D+4c_f\right)\hfill & \hfill -\mu l{c}_A\hfill & \hfill \mu l{c}_B\hfill & \hfill {\displaystyle \bar{r}}\left(c_C+2c_f\right)\hfill & \hfill -{\displaystyle \bar{r}}\left(c_D+2c_f\right)\hfill \\ \hfill {\displaystyle \bar{r}}(c_A+2c_f)\hfill & \hfill 0\hfill & \hfill c_A+2c_f\hfill & \hfill 0\hfill & \hfill -c_f\hfill & \hfill -c_f\hfill \\ \hfill -{\displaystyle \bar{r}}(c_B+2c_f)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill c_B+2c_f\hfill & \hfill -c_f\hfill & \hfill -c_f\hfill \\ \hfill 0\hfill & \hfill {\displaystyle \bar{r}}(c_C+2c_f)\hfill & \hfill -c_f\hfill & \hfill -c_f\hfill & \hfill c_C+2c_f\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -{\displaystyle \bar{r}}(c_D+2c_f)\hfill & \hfill -c_f\hfill & \hfill -c_f\hfill & \hfill 0\hfill & \hfill c_D+2c_f\hfill \end{array}\right]\text{.}$$(8)

The nominal values of the parameters, taken from [4], are: k _p  = 16 × 10⁶ kg·s⁻²; k _f   =  7 × 10⁶ kg·s⁻²; γ ₁  =  0.9; γ ₂ = 0.8; r ₁ = 75 mm; r ₂ = 120 mm; l  = 12.5 mm; c _A  = c _B  = c _C  = c _D  = 4 kg·s⁻¹; c _f   = 0.1 kg·s⁻¹.

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 is investigated by calculating the eigenvalues λ of the equations of motion (4) solving the following characteristic equation$$\det \left({\lambda }^2\text{M}+\lambda \text{C}+\text{K}\right)=0\text{.}$$(9)

Because the system (4) is a 6 DOF second-order 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):

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

2. ${\text{U}}_0$ is unstable if $\exists \mathit{i}\in (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 self-excited 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 (ξ ₁,…, ξ _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},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 α ={ α ₁,…, α _r  } a multi-index 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 non-intrusive 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 non-intrusive spectral projection (NISP) and the regression methods are the most used non-intrusive 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$${\displaystyle \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 ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\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.

4.2 Multi-element generalized polynomial chaos

This approach consists in breaking down the random space into m non-intersecting 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},i=1,\dots ,r;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 ME-gPC 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 ME-gPC, 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 low-rank 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 Low-rank 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{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 low-rank 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 m-norm 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 m-norm, the sequence of nested sets ${\mathcal{A}}_{\mathit{m}}^{\mathit{r},\mathit{p}}(\mathit{p}\in \mathbb{N})$ always converges to the set ${\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},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{{\max }_{1\le j\le r}{S}_j^T-S_i^T}{{\displaystyle \sum _{k=1}^r}{S}_k^T},i=\mathrm{1,\dots },r\text{.}$$(30)

In the above equation, $S_i^T$ is the PC-based 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 non-zero 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 ² defined as$$R^2=1-\frac{{\text{Err}}_{\text{emp}}}{\widehat{{\displaystyle \mathit{V}}}[X]}\text{,}$$(38)where ${\displaystyle \widehat{\mathit{V}}}[X]$ is the variance of X (ξ ^(q)):$${\displaystyle \widehat{\mathit{V}}}[X]=\frac{1}{Q-1}{\displaystyle \sum _{q=1}^Q}{(X({\xi }^{(q)})-\overline{\mathit{X}})}^2\text{with}{\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 leave-one-out error [36], which is based on a sum of squared predicted residuals Δ⁽ⁱ⁾ defined hereafter, may be useful to avoid this drawback. A predicted residual expresses the difference between the deterministic evaluation X(ξ ⁽ⁱ⁾) 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 (ξ ⁽¹⁾,…, ξ ^(Q))  ξ ⁽ⁱ⁾ (that is the original NED from which the observation ξ ⁽ⁱ⁾ 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 leave-one-out 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 leave-one-out 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 ² equivalent to that of the empirical error, R ², may be defined for the leave-on-out error:$$S^2=1-\frac{{\text{Err}}_{\text{LOO}}}{\widehat{\mathit{V}}[X]}\text{.}$$(43)

The two coefficients R ² and S ² 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

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)).

2. 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 low-rank index sets (or the coefficient m used for the m-norm of truncation for isotropic and anisotropic hyperbolic index sets), as well as the target accuracy $S_{\text{target}}^2$ and two thresholds ϵ ₁ and ϵ ₂.

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}\mathbf{s}\mathbf{t}\mathbf{e}\mathbf{p}$  (addition step):

1. With low-rank index sets: set the interaction order value to $j\mathrm{}=\min (\mathrm{}p\mathrm{},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}\}$ , 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 ².

   With isotropic (resp. anisotropic) index sets: for each term from the candidate set ${\mathcal{C}}_m^p=\{\alpha \in {\mathbb{N}}^r:p-1\le |\left|\alpha \right|{|}_m\le p\}$ (resp. ${\mathcal{C}}_{m,w}^p=\{\alpha \in {\mathbb{N}}^r:p-1\le |\left|\alpha \right|{|}_{m,w}\le p\}$), add it to the set ${\mathcal{A}}_m^{p-1}$ (resp. ${\mathcal{A}}_{m,w}^{p-1}$) 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 ², that is$$\mathrm{\Delta }{R}^2=R^2-R_0^2\ge {ϵ}_1\text{.}$$(44)

If inequality (44) is not respected the candidate is discarded.

With low-rank 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):

1. 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 ² in each case.

2. 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 ², i.e. if the suppression of the candidate leads to$$\mathrm{\Delta }{R}^2=R_0^2-R^2<{ϵ}_2\text{.}$$(45)

If inequality (45) is not respected the candidate is kept.

Let ${\mathcal{A}}^p$  be the final truncation set for low-rank 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 low-rank 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 low-rank 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 leave-one-out error $S_0^2$ is larger than the target value $S_{\mathrm{t}\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, ME-gPC 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 , γ ₁, γ ₂, r ₁, r ₂ 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.

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 non-zero eigenvalues (modes 1, 2, 3 and 4 with complex conjugates). Analyzing the eigenvectors corresponding to the 4 quasi-zero eigenvalues we could see that these modes are mostly constituted of wobbling and nodal-diameter 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 well-known 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 λ ₂ 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 ME-gPC 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 ME-gPC), ${\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, ME-gPC, and the three sparse PC methods that are respectively referred to as S-gPC-LR (sparse PC with low-rank index sets), S-gPC-IH (sparse PC with isotropic hyperbolic index sets) and S-gPC-AH (sparse PC with anisotropic hyperbolic index sets).

The results obtained with gPC and ME-gPC 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 ME-gPC 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 ME-gPC 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 S-gPC-LR (with low-rank 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 m-norm have been used: in the isotropic case (S-gPC-IH), the value of m was set to 0.6, a higher value leading to an excessive number of CDM; in the anisotropic case (S-gPC-AH), 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 ME-gPC, 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 S-gPC-LR and S-gPC-AH are still lower than 5%, and lower than 7% for S-gPC-IH. For r = 8, the efficiency of the sparse methods is particularly visible as the number of CDM exceeds the limit of 10 000 with ME-gPC, while the relative error becomes excessive with gPC. With the S-gPC-LR, 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 S-gPC-IH and S-gPC-AH 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 S-gPC-LR (<5% with S-gPC-AH and <7% with S-gPC-AH). 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 S-gPC-LR method whatever the number of uncertain parameters. However, the computational cost increases substantially with this method compared to S-gPC-IH and S-gPC-AH. Between the two sparse PC methods with hyperbolic index sets, i.e. S-gPC-IH and S-gPC-AH, the accuracy is seen significantly improved with the S-gPC-AH 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 S-gPC-AH 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 S-gPC-AH (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 S-gPC-AH 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 S-gPC-AH 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.

Fig. 3 Real parts of the eigenvalues of modes 1 and 2 using MCDM (black) and S-gPC-AH (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 ME-gPC.

To overcome the limitations of the classic MC approach which is often too costly for large industrial systems, gPC and ME-gPC 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 ME-gPC 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 low-rank 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 low-rank 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 low-rank 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

All Tables

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 S-gPC-AH (gray) as functions of the friction coefficient μ.
In the text

	Fig. A1 Algorithm applied to build a sparse polynomial chaos expansion with low-rank 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