Reliability-Based Robust Design of Mechanical Components with Correlated Failure Modes Based on Moment Method

Abstract

There are usually several potential failure modes in mechanical components. The conventional model for the system analysis is built under the condition that all failure modes are independent of each other. However, in engineering practice, failure modes are mostly dependent due to the fact that the elements involved in each failure mode are closely interrelated. System reliability analysis and evaluation simply conducted under independent assumption often result in excessive errors or even wrong conclusion. A novel method to evaluate system reliability of mechanical components with multiple failure modes based on moment method are proposed here. Firstly, the moment-based reliability and reliability sensitivity analysis method is proposed with independent assumption. Secondly, the proposed method is deduced by taking the correlation between failure modes into account and the correlation model is established with the copula function, which is proved to be a useful tool to model nonlinear correlation with marginal distributions. The robust design is performed as a biobjection optimization process based on the reliability sensitivity. The numerical examples show that the applied procedure is able to efficiently consider various failure modes of mechanical components in probabilistic assessment and reliability-based optimization.

1. Introduction

Structural engineers aim to better design mechanical components by reducing the cost and by increasing the performance, without compromising the structural safety. Standard of reliability has become one of the important indicators to evaluate the product quality. The aim of computational reliability analysis is to assess the reliability of a mechanical structure by considering how its performance is affected by random variations and uncertainties. The objective is to compute an estimate of the probability that the performance of the structure will satisfy a given requirement [1, 2]. Moreover, there are usually multiple potential failure modes in which failure can occur. These failure modes are often correlated because they depend on the same uncertain variables. The accuracy of the estimations highly depends on the ability to model the joint failure surface, which is established to represent the correlation between failure modes [3–5]. Thus, an efficient method to evaluate the reliability of structures under multiple failure modes is required for practical engineering design.

Research on the theory of this field has been carried out in recent decades. Cornell proposed the bounds theory on the reliability of structural systems which could be applied to obtain the rough estimation of the system reliability [6]. Ditlevsen developed the bound theory and proposed the narrow reliability bounds [7]. However, the bound theory could only obtain the reliability intervals but not the determined value. Yu et al. established the system reliability model by considering the primary and secondary failure modes of the mechanical components [8]. System reliability is fully and accurately expressed by one-order probabilities and two- and three-order joint probabilities in [9]. Probabilistic network evaluation technique is proposed to solve the system reliaiblty problem of large structure system [10]. Zhao and Ang proposed a method based on moment approximations for structural system reliability assessment which is applicable to both series and nonseries systems [11]. The work in [12] used a matrix-based system reliability (MSR) method to estimate the probabilities of complex system events by simple matrix calculations.

In reality, most of the mechanical components are affected by uncertainties, which may be caused by inherent variability (e.g., nonrepeatability of the manufacturing process, defects, etc.) or by a lack of information on the model (e.g., uncertainties on loading, deterioration processes, etc.) [13]. Uncertain variables existing in the mechanical structures may have different influence on the reliability to varying degrees. During the optimization of mechanical components, the variables that have most of the influence on the reliability of the structure have received significant attention [14–17]. Reliability sensitivity analysis has been used to obtain the change rate of a structure response due to the random inputs, which is usually solved by partial derivative or numerical method. It can rank the distribution parameters of the design variables and guide the reliability design optimization. Generally, multiple failure modes may occur when mechanical components fail. Thus, it is essential to take into account not only the potential failure modes, but also the correlation of these failure modes. So far, only a few studies have been conducted to obtain the sensitivity index for a model with correlated input. As a tool to establish a joint distribution function from its marginal distributions, copula functions are often adopted to study the correlation problems [18, 19]. Copulas are useful tools for modeling dependence among the components. They provide a way of specifying joint distributions if only the marginal distributions are known. In terms of reliability problem with multiple failure modes, we can obtain a multivariate distribution for modeling joint behavior of failure modes using the marginal distributions of each failure mode and the copula function.

Structural optimization aims at determining an optimal design according to certain criteria while meeting specific constraints [20]. Reliability-based optimization provides a well-suited framework to determine the optimal design considering uncertainties. In reliability-based robust design, the robustness of the structure reliability is searched during the optimization process to minimize the variability of output performance function when subjected to uncertainties in design parameters.

This paper presents a novel reliability optimization method based on moment method for mechanical components with multiple failure modes. The reliability analysis of failure modes of mechanical component is performed with moment method. The copula function is used to model the nonlinear correlation between failure modes. Then, the system reliability of multiple failure modes is analyzed according to the theory of probability with the copula-based joint probability distribution function. The numerical models of the reliability sensitivity with respect to the mean and the variance of the material design parameters are established. The reliability-based optimization is finally conducted by minimizing the reliability sensitivity of the design variables. The proposed method allows considering the nonlinear correlation between failure modes of mechanical component, which improves the precision of the system reliability problem.

Section 2 provides the reliability-based sensitivity analysis considering independent failure modes based on the moment method. Section 3 reformulates the method for reliability analysis at the system level with copula functions. Section 4 introduces the reliability-based optimization process by considering both the probabilistic constraints and the reliability sensitivity as the subobjective. Section 5 applies the newly derived method to three system problems with multiple failure modes comparing its accuracy to sampling methods. Section 6 provides some conclusions on the proposed method.

2. Reliability Analysis considering Independent Failure Modes

The assumption of independent failure modes could reduce the difficulty of the problem and it gives the estimation of the reliability of mechanical components with multiple failure modes in an efficient way, as well as the reliability sensitivity of the random parameters with respect to the structural reliability. The reliability of a mechanical component with multiple failure modes could be expressed as

$R=P\left(G>\mathrm{0}\right)={{\displaystyle \int }}_{\mathrm{0}}^{\infty }{f}_G\left(G\right)\text{d}G,$

(1)

where $\mathbf{X}={[X_{\mathrm{1}},X_{\mathrm{2}},\dots ,X_n]}^{\text{T}}$ denotes the random variable vector and f _G (G) defines the joint probability density function.

According to the reliability study of different forms, the failure modes could be expressed in terms of several limit state functions as a set of equations as follows:

$\begin{array}{c}{g}_{\mathrm{1}}\left(\mathbf{X}\right)=\mathrm{0}\\ g_{\mathrm{2}}\left(\mathbf{X}\right)=\mathrm{0}\\ ⋮\\ g_m\left(\mathbf{X}\right)=\mathrm{0},\end{array}$

(2)

where g _i (X) represents the performance function of the ith failure mode.

According to the reliability theory, the reliability index of each failure mode is defined as

$\begin{array}{c}{\beta }_{\mathrm{2}Mi}=\frac{\text{E}\left[g_i\left(\mathbf{X}\right)\right]}{\sqrt{\mathrm{Var}\left[g_i\left(\mathbf{X}\right)\right]}}=\frac{{\mu }_{g_i}}{{\sigma }_{g_i}},\\ R_i=\mathrm{\Phi }\left({\beta }_{\mathrm{2}Mi}\right),\end{array}$

(3)

where μ _{g i} , σ _{g i} , and β_2Mi are the mean, the variance, and the reliability index, respectively. Ф(·) is the standard normal distribution function.

In the case that the random variables are normally distributed, (3) could accurately estimate the reliability. However, the distribution types of the random variables are unknown on most occasions and the complete probabilistic information needed is unable to obtain. To overcome this problem, the moment method could estimate the reliability index by making use of the first few moments of random variables. One of the moment methods called high-order moment standardization technique is adopted here and the reliability index of each failure mode could be expressed as follows according to the theory in [11]:

${\beta }_{\mathrm{4}Mi}=\frac{\mathrm{3}\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right){\beta }_{\mathrm{2}Mi}+{\alpha }_{\mathrm{3}gi}\left({\beta }_{\mathrm{2}Mi}^{\mathrm{2}}-\mathrm{1}\right)}{\sqrt{\left(\mathrm{9}{\alpha }_{\mathrm{4}gi}-\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}\right)\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)}},$

(4)

where ${\alpha }_{\mathrm{3}gi}={\theta }_{gi}/{\sigma }_{gi}^{\mathrm{3}}$ , ${\alpha }_{\mathrm{4}gi}={\eta }_{gi}/{\sigma }_{gi}^{\mathrm{4}}$ represent the coefficients of skewness and kurtosis of the ith failure mode, respectively. θ _gi and η _gi are the third and the fourth central moments of the ith failure mode. According to the perturbation technique, θ _gi and η _gi have the following expressions:

$\begin{array}{c}{\theta }_{gi}=\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\mu }}_{\mathrm{4}}\left(\mathbf{X}\right){\left[\frac{\partial g\left(\mathbf{X}\right)}{\partial {\mathbf{X}}^{\text{T}}}\right]}^{\left[\mathrm{2}\right]}\frac{{\partial }^{\mathrm{2}}g\left(\mathbf{X}\right)}{\partial {\mathbf{X}}^{\text{T}\mathrm{2}}},\\ {\eta }_{gi}={\mathit{\mu }}_{\mathrm{4}}\left(\mathbf{X}\right){\left[\frac{\partial g\left(\mathbf{X}\right)}{\partial {\mathbf{X}}^{\text{T}}}\right]}^{\left[\mathrm{4}\right]},\end{array}$

(5)

where $\mathrm{0}<\left|\epsilon \right|\ll \mathrm{1}$ is a small parameter and ${\mathit{\mu }}_{\mathrm{4}}(\mathbf{X})$ denotes the fourth moment vector.

The reliability of each failure mode can be obtained with the reliability index β_4Mi:

$R_i=\mathrm{\Phi }\left({\beta }_{\mathrm{4}Mi}\right).$

(6)

Generally, the reliability of a mechanical component with m correlated failure modes could be expressed as

$\begin{array}{c}R=P\left\{g_{\mathrm{1}}\left(\mathbf{X}\right)>\mathrm{0}\cap g_{\mathrm{2}}\left(\mathbf{X}\right)>\mathrm{0}\cap \cdots g_n\left(\mathbf{X}\right)>\mathrm{0}\right\}\\ ={{\displaystyle \int }}_{\mathrm{0}}^{\infty }\cdots {{\displaystyle \int }}_{\mathrm{0}}^{\infty }{f}_G\left(g_{\mathrm{1}},g_{\mathrm{2}},\dots ,g_m\right)\text{d}{g}_{\mathrm{1}}\text{d}{g}_{\mathrm{2}}\cdots \text{d}{g}_n,\end{array}$

(7)

where $f_G(g_{\mathrm{1}},g_{\mathrm{2}},\dots ,g_m)$ is the joint probability density function of random variables of each failure mode.

Unlike large-scale structural system, mechanical components usually fail with less failure modes. In this case, the failure of the mechanical component with multiple failure modes could be considered as a series system reliability problem under the independence assumption, the reliability of which is as follows:

$R={\displaystyle \prod }_{i=\mathrm{1}}^m{R}_i.$

(8)

3. The Reliability-Based Sensitivity Analysis with Independent Failure Modes

The reliability-based sensitivity refers to the partial derivative of the structural reliability with respect to random variables. Based on (4), the reliability sensitivity could be deduced by direct differential method. Considering the arbitrary distributed random variables, the reliability sensitivity with respect to the mean and the variance of random variables can be obtained as follows:

$\begin{array}{c}\frac{\text{D}R}{\text{D}{\stackrel{-}{\mathbf{X}}}^{\text{T}}}={\displaystyle \sum }_{i=\mathrm{1}}^n\{\frac{\partial R_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\mu }_{gi}}\frac{\partial {\mu }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\\ \hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{1em}+\left[\frac{\partial R_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}\left(\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}+\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\sigma }_{gi}}\right)\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}\times \frac{\partial {\sigma }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\phantom{\frac{\partial R_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\mu }_{gi}}\frac{\partial {\mu }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}}\}\hspace{1em}\times \left(\frac{{{\displaystyle \prod }}_{j=\mathrm{1}}^n{R}_j\left({\beta }_{\mathrm{4}Mj}\right)}{R_i\left({\beta }_{\mathrm{4}Mi}\right)}\right),\\ \frac{\text{D}R}{\text{DVar}\left(\mathbf{X}\right)}={\displaystyle \sum }_{i=\mathrm{1}}^n\frac{\partial R_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}\left(\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}+\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\sigma }_{gi}}\right)\hspace{1em}\hspace{1em}\times \frac{\partial {\sigma }_{gi}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}\left(\frac{{{\displaystyle \prod }}_{j=\mathrm{1}}^n{R}_j\left({\beta }_{\mathrm{4}Mj}\right)}{R_i\left({\beta }_{\mathrm{4}Mi}\right)}\right),\end{array}$

(9)

where

$\begin{array}{c}\frac{\partial R_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}=\phi \left({\beta }_{\mathrm{4}Mi}\right),\\ \frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}=\frac{\mathrm{3}\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)+\mathrm{2}{\alpha }_{\mathrm{3}gi}{\beta }_{\mathrm{2}Mi}}{\sqrt{\left(\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}{\alpha }_{\mathrm{4}gi}+\mathrm{9}\right)\left(\mathrm{1}-{\alpha }_{\mathrm{4}gi}\right)}},\\ \frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\mu }_{gi}}=\frac{\mathrm{1}}{{\sigma }_{gi}},\\ \frac{\partial {\mu }_{gi}}{\partial {\mathbf{X}}^{\text{T}}}=\left[\frac{\partial g_i\left(\stackrel{-}{\mathbf{X}}\right)}{\partial X_{\mathrm{1}}}\frac{\partial g_i\left(\stackrel{-}{\mathbf{X}}\right)}{\partial X_{\mathrm{2}}}\cdots \frac{\partial g_i\left(\stackrel{-}{\mathbf{X}}\right)}{\partial X_n}\right],\\ \frac{\partial {\sigma }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\mathrm{1}}{\mathrm{2}{\sigma }_{gi}}[\frac{{\partial }^{\mathrm{2}}{g}_i}{\partial {\left({\mathbf{X}}^{\text{T}}\right)}^{\mathrm{2}}}\otimes \frac{\partial g_i}{\partial {\mathbf{X}}^{\text{T}}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left(\frac{{\partial }^{\mathrm{2}}{g}_i}{\partial {\left({\mathbf{X}}^{\text{T}}\right)}^{\mathrm{2}}}\otimes \frac{\partial g_i}{\partial {\mathbf{X}}^{\text{T}}}\right)\left({\mathbf{I}}_n\otimes {\mathbf{U}}_{n^{\mathrm{2}}\times n^{\mathrm{2}}}\right)]\hspace{1em}\times \left({\mathbf{I}}_n\otimes \mathrm{Var}\left(\mathbf{X}\right)\right),\\ \frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}=-\frac{{\mu }_{gi}}{{\sigma }_{gi}^{\mathrm{2}}},\\ \frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\sigma }_{gi}}\hspace{1em}=\frac{\left({\alpha }_{\mathrm{3}gi}\left(\mathrm{3}-\mathrm{5}{\beta }_{\mathrm{2}Mi}^{\mathrm{2}}\right)-\mathrm{15}{\alpha }_{\mathrm{4}gi}{\beta }_{\mathrm{2}Mi}+\mathrm{3}{\beta }_{\mathrm{2}Mi}\right)/{\sigma }_{gi}}{\sqrt{\left(\mathrm{9}{\alpha }_{\mathrm{4}gi}-\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}\right)\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)}}\hspace{1em}\hspace{1em}-\frac{\mathrm{1}}{\mathrm{2}}\frac{\left[\mathrm{3}{\beta }_{\mathrm{2}Mi}\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)+{\alpha }_{\mathrm{3}gi}\left({\beta }_{\mathrm{2}Mi}^{\mathrm{2}}-\mathrm{1}\right)\right]}{\sqrt{{\left(\mathrm{9}{\alpha }_{\mathrm{4}gi}-\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}\right)}^{\mathrm{3}}{\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)}^{\mathrm{3}}}}\hspace{1em}\hspace{1em}\times \frac{\left(\mathrm{50}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}{\alpha }_{\mathrm{4}gi}-\mathrm{30}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{72}{\alpha }_{\mathrm{4}gi}^{\mathrm{2}}+\mathrm{72}{\alpha }_{\mathrm{4}gi}\right)}{{\sigma }_{gi}},\\ \frac{\partial {\sigma }_{gi}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\mathrm{1}}{\mathrm{2}{\sigma }_{gi}}\left[\frac{\partial {\stackrel{-}{g}}_i}{\partial \mathbf{X}}\otimes \frac{\partial {\stackrel{-}{g}}_i}{\partial \mathbf{X}}\right].\end{array}$

(10)

Substituting (10) into (9), the reliability sensitivity of random variables $\text{D}R/\text{D}{\stackrel{-}{\mathbf{X}}}^{\text{T}}$ and $\text{D}R/\text{DVar}(\mathbf{X})$ can be obtained.

Due to the fact that different kinds of random parameters are involved in the design process, such as material design parameter and load parameters, it is necessary to transform these parameters into a set of dimensionless parameters. Thus, the reliability sensitivity of random parameters could be able to be compared and the most sensitive parameters could be recognized.

The nondimensional method of the reliability sensitivity with respect to the mean and the variance can be expressed as follows:

$\begin{array}{c}{\alpha }_i=\frac{\text{D}{R}_f}{\text{D}{\stackrel{-}{X}}_i}\frac{{\sigma }_X}{P_f},\\ {\eta }_i=\frac{\text{D}{R}_f}{\text{D}{\sigma }_{X_i}}\frac{{\mathrm{Var}}_X}{P_f},\end{array}$

(11)

where σ _X and Var _X represent the standard deviation and the variance of the random variables, respectively, and P _f is the failure probability.

4. The Reliability-Based Sensitivity Analysis with Correlated Failure Modes

In this section, the correlation of different failure modes is taken into consideration. Actually, because the degree of correlation of failure modes is different, the failure probability of one failure mode may be influenced by another. When two failure modes are highly correlated, the failure of one failure mode would lead to another failure. The two failure modes could be considered as independent modes if they are slightly correlated. From this point of view, it is necessary to take into account the correlation form of failure modes and also the correlation degree.

4.1. Copula Definition

The copula of a multivariate distribution describes not only the correlations of the random variables, but also the dependence structure. It is uncoupled from the marginal distributions which can be modeled as empirical distributions or fitted standard distributions as usual [21].

Let F _X (x) and F _Y (y) denote the marginal distribution functions of the variables X and Y, respectively. The joint distribution function F_X,Y(x,y) can be expressed as

$F_{X,Y}\left(x,y\right)=C\left[F_X\left(x\right),F_Y\left(y\right)\right],$

(12)

where C(u,v) is the copula function.

If F _X (x) and F _Y (y) are continuous functions, C(u,v) is unique. Otherwise, if F _X (x) and F _Y (y) are univariate functions and C(u,v) is a copula function, then F_X,Y(x,y) is a bivariate joint distribution function with marginal F _X (x) and F _Y (y).

Frees and Valdez provide a number of copulas [22]. Following are some commonly used bivariate copula functions.

Gaussian copula is

$C\left(u,v\right)=H\left[{\mathrm{\Phi }}^{-\mathrm{1}}\left(u\right),{\mathrm{\Phi }}^{-\mathrm{1}}\left(v\right)\right],\hspace{1em}-\mathrm{1}\le \alpha \le \mathrm{1},$

(13)

where Φ⁻¹(·) is the inverse of the univariate standard normal distribution and H(x,y) is a standard normal distribution function.

Gumbel copula is

$C\left(u,v\right)=\exp \left(-{\left[{\left(-\ln u\right)}^{\mathrm{1}/\alpha }+{\left(-\ln v\right)}^{\mathrm{1}/\alpha }\right]}^{\alpha }\right),\hspace{1em}\mathrm{0}<\alpha \le \mathrm{1}.$

(14)

Frank copula is

$C\left(u,v\right)=-\frac{\mathrm{1}}{\alpha }\ln \left[\mathrm{1}+\frac{\left({\text{e}}^{-\alpha u}-\mathrm{1}\right)\left({\text{e}}^{-\alpha v}-\mathrm{1}\right)}{{\text{e}}^{-\alpha }-\mathrm{1}}\right],\hspace{1em}\alpha \ne \mathrm{0}.$

(15)

For more formal definitions of copula functions, the reader is referred to [23].

Generally, the failure modes of the mechanical components, such as yield, fracture, cumulative fatigue, corrosion, creep, wear, and deformation exceeding limits, are basically positive correlated; that is, one failure mode would lead to an increased failure probability of another one [24]. Besides, the first derivative of the copula function with respect to random variables is often needed during the reliability optimization process. However, many copula functions are expressed in a complex form that leads to a more complex expression of the first derivative (or even unavailable). Thus, the selection of the proper copula function describing the correlation between failure modes becomes a considerable problem. As one of the Archimedes copula functions, the Gumbel copula is proved to be able to model the positive correlation between random variables, the expression of which is also concise. In this paper, the Gumbel copula is selected to model the joint probability distribution function of the possible failure modes of mechanical components, and the reliability optimization problem could then be conducted.

The expression of the first derivative of the Gumbel copula is as follows:

$\begin{array}{c}\frac{\partial C_G}{\partial u}=-(\phantom{{\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}}{\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }\hspace{1em}\hspace{1em}\hspace{1em}\times {\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }})\hspace{1em}\times {\left(u\ln \left(u\right)\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)\right)}^{-1},\\ \frac{\partial C_G}{\partial v}=-(\phantom{{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}}{\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }\hspace{1em}\hspace{1em}\times {\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }})\times {\left(v\ln \left(v\right)\left({\left(-\ln \left(u\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(v\right)\right)}^{\mathrm{1}/\alpha }\right)\right)}^{-1}.\end{array}$

(16)

The most commonly used copula functions are bivariant copula functions. Especially, due to the characteristics of the Archimedes copula functions, the bivariant Archimedes copula functions are more widely used to model correlations. In fact, any N-dimension Archimedes copula function could be deduced from a bivariant one [23]. Thus, the reliability problem with multiple failure modes could be transformed into the bivariant form

$\begin{array}{c}C\left(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}}\right)=C\left(C\left(u_{\mathrm{1}},u_{\mathrm{2}}\right),u_{\mathrm{3}}\right),\\ C\left(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}},u_{\mathrm{4}}\right)=C\left(C\left(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}}\right),u_{\mathrm{4}}\right),\\ C\left(u_{\mathrm{1}},u_{\mathrm{2}},\dots ,u_{N-\mathrm{1}},u_N\right)=C\left(C\left(u_{\mathrm{1}},u_{\mathrm{2}},\dots ,u_{N-\mathrm{1}}\right),u_N\right).\end{array}$

(17)

Firstly, the bivariant copula function is adopted when considering only two failure modes for the fact that the theory and formula of which are relatively simple. In this case, the joint probability distribution function of the two failure modes could be directly built with the selected bivariant copula function. The reliability sensitivity of the random variables with respect to the mean and the variance could then be obtained by the matrix differential technology. Then, the method is extended to the case of N correlated failure modes.

4.2. The Case of Two Correlated Failure Modes

According to the theory of probability, the failure probability of a mechanical component with two failure modes could be expressed as

$\begin{array}{c}P=P\left(g_{\mathrm{1}}\left(\mathbf{X}\right)\le \mathrm{0}\cup g_{\mathrm{2}}\left(\mathbf{X}\right)\le \mathrm{0}\right)\\ =P\left(g_{\mathrm{1}}\left(\mathbf{X}\right)\le \mathrm{0}\right)+P\left(g_{\mathrm{2}}\left(\mathbf{X}\right)\le \mathrm{0}\right)\\ \hspace{1em}-P\left(g_{\mathrm{1}}\left(\mathbf{X}\right)\le \mathrm{0}\cap g_{\mathrm{2}}\left(\mathbf{X}\right)\le \mathrm{0}\right)\\ =P_{\mathrm{1}}+P_{\mathrm{2}}-C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),\end{array}$

(18)

where $\mathbf{X}={[X_{\mathrm{1}},X_{\mathrm{2}},\dots ,X_n]}^{\text{T}}$ is the random variable vector and the reliability sensitivity of the random variables with respect to the mean and the variance is as follows:

$\begin{array}{c}\frac{\partial P}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\partial P_{\mathrm{1}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}+\frac{\partial P_{\mathrm{2}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}-\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}},\\ \frac{\partial P}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\partial P_{\mathrm{1}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}+\frac{\partial P_{\mathrm{2}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}-\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}.\end{array}$

(19)

Based on the first four-moment and the high-order moment standardization technique, the reliability sensitivity is obtained as

$\begin{array}{c}\frac{\partial P_i}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\partial P_i\left({\beta }_{4Mi}\right)}{\partial {\beta }_{4Mi}}\frac{\partial {\beta }_{4Mi}}{\partial {\beta }_{2Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\mu }_{gi}}\frac{\partial {\mu }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\\ \hspace{1em}+\left[\frac{\partial P_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}\left(\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}+\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\sigma }_{gi}}\right)\right]\frac{\partial {\sigma }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}},\end{array}$

(20)

$\begin{array}{c}\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\left(\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{1}}}\frac{\partial P_{\mathrm{1}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}+\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{2}}}\frac{\partial P_{\mathrm{2}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\right),\end{array}$

(21)

$\begin{array}{c}\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{1}}}=-(\phantom{{\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}}{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }\times {\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }\times {\text{e}}^{-{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }})\hspace{1em}\times (\phantom{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}{P}_{\mathrm{1}}\ln \left(P_{\mathrm{1}}\right){\hspace{0.5em}\times \left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right))}^{-1},\end{array}$

(22)

$\begin{array}{c}\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{2}}}=-(\phantom{{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}}{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }\times {\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }})\times (\phantom{\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right)}{P}_{\mathrm{2}}\ln \left(P_{\mathrm{2}}\right)\times {\left({\left(-\ln \left(P_{\mathrm{1}}\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_{\mathrm{2}}\right)\right)}^{\mathrm{1}/\alpha }\right))}^{-1},\end{array}$

(23)

$\frac{\partial P_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}=-\phi \left(-{\beta }_{\mathrm{4}Mi}\right),$

(24)

$\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}=\frac{\mathrm{3}\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)+\mathrm{2}{\alpha }_{\mathrm{3}gi}{\beta }_{\mathrm{2}Mi}}{\sqrt{\left(\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}{\alpha }_{\mathrm{4}gi}+\mathrm{9}\right)\left(\mathrm{1}-{\alpha }_{\mathrm{4}gi}\right)}},$

(25)

$\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\mu }_{gi}}=\frac{\mathrm{1}}{{\sigma }_{gi}},$

(26)

$\frac{\partial {\mu }_{gi}}{\partial {\mathbf{X}}^{\text{T}}}=\left[\frac{\partial g_i\left(\stackrel{-}{\mathbf{X}}\right)}{\partial X_{\mathrm{1}}}\frac{\partial g_i\left(\stackrel{-}{\mathbf{X}}\right)}{\partial X_{\mathrm{2}}}\right],$

(27)

$\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}=-\frac{{\mu }_{gi}}{{\sigma }_{gi}^{\mathrm{2}}},$

(28)

$\begin{array}{c}\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\sigma }_{gi}}=\frac{\left({\alpha }_{\mathrm{3}gi}\left(\mathrm{3}-\mathrm{5}{\beta }_{\mathrm{2}Mi}^{\mathrm{2}}\right)-\mathrm{15}{\alpha }_{\mathrm{4}gi}{\beta }_{\mathrm{2}Mi}+\mathrm{3}{\beta }_{\mathrm{2}Mi}\right)/{\sigma }_{gi}}{\sqrt{\left(\mathrm{9}{\alpha }_{\mathrm{4}gi}-\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}\right)\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)}}\\ \hspace{1em}-\frac{\mathrm{1}}{\mathrm{2}}\frac{\left[\mathrm{3}{\beta }_{\mathrm{2}Mi}\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)+{\alpha }_{\mathrm{3}gi}\left({\beta }_{\mathrm{2}Mi}^{\mathrm{2}}-\mathrm{1}\right)\right]}{\sqrt{{\left(\mathrm{9}{\alpha }_{\mathrm{4}gi}-\mathrm{5}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{9}\right)}^{\mathrm{3}}{\left({\alpha }_{\mathrm{4}gi}-\mathrm{1}\right)}^{\mathrm{3}}}}\hspace{1em}\\ \times \frac{\left(\mathrm{50}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}{\alpha }_{\mathrm{4}gi}-\mathrm{30}{\alpha }_{\mathrm{3}gi}^{\mathrm{2}}-\mathrm{72}{\alpha }_{\mathrm{4}gi}^{\mathrm{2}}+\mathrm{72}{\alpha }_{\mathrm{4}gi}\right)}{{\sigma }_{gi}},\end{array}$

(29)

$\frac{\partial {\sigma }_{gi}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\mathrm{1}}{\mathrm{2}{\sigma }_{gi}}[\frac{{\partial }^{\mathrm{2}}{g}_i}{\partial {\left({\mathbf{X}}^{\text{T}}\right)}^{\mathrm{2}}}\otimes \frac{\partial g_i}{\partial {\mathbf{X}}^{\text{T}}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}+\left(\frac{{\partial }^{\mathrm{2}}{g}_i}{\partial {\left({\mathbf{X}}^{\text{T}}\right)}^{\mathrm{2}}}\otimes \frac{\partial g_i}{\partial {\mathbf{X}}^{\text{T}}}\right)\left({\mathbf{I}}_n\otimes {\mathbf{U}}_{n^{\mathrm{2}}\times n^{\mathrm{2}}}\right)]\hspace{1em}\times \left({\mathbf{I}}_n\otimes \mathrm{Var}\left(\mathbf{X}\right)\right),$

(30)

$\frac{\text{D}{P}_i}{\text{DVar}\left(\mathbf{X}\right)}=\frac{\partial P_i\left({\beta }_{\mathrm{4}Mi}\right)}{\partial {\beta }_{\mathrm{4}Mi}}\hspace{1em}\times \left(\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\beta }_{\mathrm{2}Mi}}\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}+\frac{\partial {\beta }_{\mathrm{4}Mi}}{\partial {\sigma }_{gi}}\right)\frac{\partial {\sigma }_{gi}}{\partial \mathrm{Var}\left(\mathbf{X}\right)},$

(31)

$\begin{array}{c}\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{1}}}\frac{\partial P_{\mathrm{1}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}\\ \hspace{1em}+\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{2}}}\frac{\partial P_{\mathrm{2}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)},\end{array}$

(32)

$\frac{\partial {\beta }_{\mathrm{2}Mi}}{\partial {\sigma }_{gi}}=-\frac{{\mu }_{gi}}{{\sigma }_{gi}^{\mathrm{2}}},$

(33)

$\frac{\partial {\sigma }_{gi}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\mathrm{1}}{\mathrm{2}{\sigma }_{gi}}\left[\frac{\partial {\stackrel{-}{g}}_i}{\partial \mathbf{X}}\otimes \frac{\partial {\stackrel{-}{g}}_i}{\partial \mathbf{X}}\right],$

(34)

where C(·) is the Gumbel copula function.

4.3. The Case of N Correlated Failure Modes

Based on the equations obtained above, the reliability analysis with N failure modes could be carried out similarly. According to the probability theory, the failure probability of a mechanical component with N failure modes could be expressed as

$\begin{array}{c}P=P\left({\displaystyle \bigcup }_{i=\mathrm{1}}^N{g}_i\left(X_{\mathrm{1}},X_{\mathrm{2}},\dots ,X_n\right)\le \mathrm{0}\right)\\ ={\displaystyle \sum }_{i=\mathrm{1}}^N{P}_i-{\displaystyle \sum }_{\mathrm{1}\le i<h\le N}C\left(P_i,P_h\right)\\ \hspace{1em}+{\displaystyle \sum }_{\mathrm{1}\le i<h<t\le N}C\left(P_i,P_h,P_t\right)-\cdots +{\left(-\mathrm{1}\right)}^{N-\mathrm{1}}\\ \hspace{1em}\times C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_N\right),\end{array}$

(35)

where C(·) is the copula function. As mentioned above, the expressions of N-dimension copula function are difficult to obtain and the parameter estimation is also a challenge. The Archimedes copulas with the characteristic that an N-dimension copula function could be represented by the bivariant copula function facilitate the correlation modeling of multiple failure modes. Therefore, the reliability sensitivity considering correlated failure modes could be expressed as follows:

$\begin{array}{c}\frac{\partial P}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}={\displaystyle \sum }_{i=\mathrm{1}}^m\frac{\partial P_i}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}-{\displaystyle \sum }_{\mathrm{1}\le i\le h\le m}\frac{\partial C\left(P_i,P_h\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\hspace{1em}+{\displaystyle \sum }_{\mathrm{1}\le i\le h\le t\le m}\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\hspace{1em}-\cdots +{\left(-\mathrm{1}\right)}^{m-\mathrm{1}}\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}},\\ \frac{\partial P}{\partial \mathrm{Var}\left(\mathbf{X}\right)}={\displaystyle \sum }_{i=\mathrm{1}}^m\frac{\partial P_i}{\partial \mathrm{Var}\left(\mathbf{X}\right)}-{\displaystyle \sum }_{\mathrm{1}\le i\le h\le m}\frac{\partial C\left(P_i,P_h\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}\hspace{1em}+{\displaystyle \sum }_{\mathrm{1}\le i\le h\le t\le m}\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}\hspace{1em}-\cdots +{\left(-\mathrm{1}\right)}^{m-\mathrm{1}}\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)},\end{array}$

(36)

where,

$\begin{array}{c}\frac{\partial C\left(P_i,P_h\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\partial C\left(P_i,P_h\right)}{\partial P_i}\frac{\partial P_i}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}+\frac{\partial C\left(P_i,P_h\right)}{\partial P_h}\frac{\partial P_h}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}},\\ \frac{\partial C\left(P_i,P_h\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\partial C\left(P_i,P_h\right)}{\partial P_i}\frac{\partial P_i}{\partial \mathrm{Var}\left(\mathbf{X}\right)}+\frac{\partial C\left(P_i,P_h\right)}{\partial P_h}\frac{\partial P_h}{\partial \mathrm{Var}\left(\mathbf{X}\right)},\\ \frac{\partial C\left(P_i,P_h,P_t\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial C\left(P_i,P_h\right)}\frac{\partial C\left(P_i,P_h\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}+\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial P_t}\frac{\partial P_t}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial C\left(P_i,P_h\right)}\times \left(\frac{\partial C\left(P_i,P_h\right)}{\partial P_i}\frac{\partial P_i}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}+\frac{\partial C\left(P_i,P_h\right)}{\partial P_h}\frac{\partial P_h}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\right)+\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial P_t}\frac{\partial P_t}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}},\\ \frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}=\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{m-\mathrm{1}}\right)}\times (\phantom{\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}\right),P_{m-\mathrm{1}}\right)}{\partial P_{m-\mathrm{1}}}\frac{\partial P_{m-\mathrm{1}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}}\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{m-\mathrm{2}}\right),P_{m-\mathrm{1}}\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{m-\mathrm{2}}\right)}\cdots (\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),P_{\mathrm{3}}\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}\times \left(\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{1}}}\frac{\partial P_{\mathrm{1}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}+\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{2}}}\frac{\partial P_{\mathrm{2}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}}\right)+\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),P_{\mathrm{3}}\right)}{\partial P_{\mathrm{3}}}\frac{\partial P_{\mathrm{3}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}})+\cdots +\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}\right),P_{m-\mathrm{1}}\right)}{\partial P_{m-\mathrm{1}}}\frac{\partial P_{m-\mathrm{1}}}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}})+\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{\mathrm{4}}\right),P_m\right)}{\partial P_m}\frac{\partial P_m}{\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}},\\ \frac{\partial C\left(P_i,P_h,P_t\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial C\left(P_i,P_h\right)}\frac{\partial C\left(P_i,P_h\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}+\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial P_t}\frac{\partial P_t}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial C\left(P_i,P_h\right)}\times (\frac{\partial C\left(P_i,P_h\right)}{\partial P_i}\frac{\partial P_i}{\partial \mathrm{Var}\left(\mathbf{X}\right)}\hspace{0.5em}+\frac{\partial C\left(P_i,P_h\right)}{\partial P_h}\frac{\partial P_h}{\partial \mathrm{Var}\left(\mathbf{X}\right)})+\frac{\partial C\left(P_i,P_h,P_t\right)}{\partial P_t}\frac{\partial P_t}{\partial \mathrm{Var}\left(\mathbf{X}\right)},\\ \frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial \mathrm{Var}\left(\mathbf{X}\right)}=\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{m-\mathrm{1}}\right)}\times (\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{m-\mathrm{2}}\right),P_{m-\mathrm{1}}\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{m-\mathrm{2}}\right)}\cdots (\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),P_{\mathrm{3}}\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}\hspace{0.5em}\times (\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{1}}}\frac{\partial P_{\mathrm{1}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)}+\frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)}{\partial P_{\mathrm{2}}}\frac{\partial P_{\mathrm{2}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)})+\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),P_{\mathrm{3}}\right)}{\partial P_{\mathrm{3}}}\frac{\partial P_{\mathrm{3}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)})+\cdots +\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}\right),P_{m-\mathrm{1}}\right)}{\partial P_{m-\mathrm{1}}}\frac{\partial P_{m-\mathrm{1}}}{\partial \mathrm{Var}\left(\mathbf{X}\right)})+\frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{m-\mathrm{1}}\right),P_m\right)}{\partial P_m}\frac{\partial P_m}{\partial \mathrm{Var}\left(\mathbf{X}\right)},\\ \frac{\partial C\left(C\left(P_i,P_h\right),P_t\right)}{\partial C\left(P_i,P_h\right)}=-{\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }\times (\phantom{\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)}C\left(P_i,P_h\right)\ln \left(C\left(P_i,P_h\right)\right){\times \left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right))}^{-1}*{\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }\times {\text{e}}^{-{\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }},\\ \frac{\partial C\left(C\left(P_i,P_h\right),P_t\right)}{\partial P_t}=-(\phantom{{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}}{\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }\hspace{0.5em}\times {\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}){\times \left(P_t\ln \left(P_t\right)\left({\left(-\ln \left(C\left(P_i,P_h\right)\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_t\right)\right)}^{\mathrm{1}/\alpha }\right)\right)}^{-1},\\ \frac{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_m\right)}{\partial C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{m-\mathrm{1}}\right)}=-\frac{{\left({\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}{C_{\Sigma }\ln \left(C_{\Sigma }\right)\left({\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }\right)}\times {\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }},\\ \frac{\partial C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{m-\mathrm{1}}\right),P_m\right)}{\partial P_m}=-\frac{{\left({\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }}{P_m\ln \left(P_m\right)\left({\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }\right)}\times {\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }{\text{e}}^{-{\left({\left(-\ln \left(C_{\Sigma }\right)\right)}^{\mathrm{1}/\alpha }+{\left(-\ln \left(P_m\right)\right)}^{\mathrm{1}/\alpha }\right)}^{\alpha }},\end{array}$

(37)

where, $\partial P_i/\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}$ is equal to (20), $\partial P_i({\beta }_{\mathrm{4}Mi})/\partial {\beta }_{\mathrm{4}Mi}$ is equal to (24), ∂β_4Mi/∂β_2Mi is equal to (25), ∂β_2Mi/∂μ _gi is equal to (26), ∂μ _gi /∂X^T is equal to (27), ∂β_2Mi/∂σ _gi is equal to (28), ∂β_4Mi/∂σ _gi is equal to (29), $\partial {\sigma }_{gi}/\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}$ is equal to (30), $\partial P_i/\partial \mathrm{Var}(\mathbf{X})$ is equal to (31), ∂β_2Mi/∂σ _gi is equal to (33), $\partial {\sigma }_{gi}/\partial \mathrm{Var}(\mathbf{X})$ is equal to (34), and $C_{\Sigma }=C(P_{\mathrm{1}},P_{\mathrm{2}},\dots ,P_{m-\mathrm{1}})$ .

According to the known probabilistic information of the random variables, the reliability sensitivity of the mechanical component with N correlated failure modes could be obtained.

4.4. Reliability Estimation with Implicit Limit State Function Problem

Structural reliability methods aim at computing the probability of failure of mechanical structures with respect to some prescribed limit state functions, which is usually resort to running an expensive-to-evaluate computational model (e.g., a finite element model) in modern engineering. The difficulty of the reliability analysis depends on the complexity of mechanical structures. For simple mechanical structures, such as a beam or an axle, it is easy to obtain the stress expression of the dangerous area with precise mechanical formula derivation. Thus it is obviously convenient to perform reliability analysis by simplifying complex structures into simple ones. These kinds of reliability problems are considered as the explicit limit state function problem; the reliability analysis and optimization could then be directly carried out. However, mechanical products normally have more complex geometric structure, the stress analysis of which can only be realized by some other techniques, such as finite element analysis. From the reliability point of view, the stress expression of the structure could not be obtained by mechanical derivation, which results in the implicit limit state function problem. In this case, because obtaining the explicit limit state functions is not straightforward, the reliability theory depending on the limit state functions can not be directly applied. Therefore, it is critical to establish the expression of the structure responses and solve the implicit limit state function problem. Among the existing analysis methods, the design of experiment (DOE) is normally used as a tool to solve the problem. The data of the input parameters and response output could be obtained by DOE, and the data are selected as samples which could be fitted into the explicit expression.

The simulation-based optimization design of the structural response includes a lot of repeated function calculations and cost significant computing resources, resulting in the fact that the computer simulation such as Monte Carlo simulation is not suitable for the optimization design problems. In order to overcome this defect, the surrogate model is widely adopted in the reliability-based robust optimization problems. Surrogate models such as quadratic response surfaces, artificial neural networks (ANN), or Kriging are then introduced as a substitute for the original model to cope with implicit limit state function problem.

Thus, the establishment of the surrogate model is a premise while conducting the reliability-based optimization design. The procedure of the modeling process is as follows.

1. Design of experiment (DOE): set the number of the samples according to the problem under study and select proper sampling method.

2. Calculate the response of each sample point with a certain calculation method. Generally, the finite element method is used to obtain the response of complex structures, noting that it is time consuming to calculate the response of each sample point when the finite element model is large.

3. Select the proper surrogate model and estimate the parameter of the model according to the original samples (input and output samples). There is generally one or several parameters in the surrogate model which needs to be accurately estimated with parameter estimation methods, such as solving linear equations system, least square method or moving least square method, and maximum likelihood approximation method.

4. Response prediction according to the established surrogate model: the response analysis process, such as finite element analysis, could be replaced by the established surrogate model.

In this paper, software, such as Pro/E, Hypermesh, MSC. Patran, MSC. Nastran, iSIGHT, and Matlab, is integrated to realize the reliability and reliability sensitivity analysis with respect to the implicit limit state function problem. Pro/E is used to establish the parametric model of the mechanical structure and the design parameters are parameterized. The finite element model is built with Hypermesh. Then, Pro/E, MSC. Patran, and MSC. Nastran are integrated by using iSIGHT. Based on the input samples generated by Latin hypercube sample technique, the DOE study could be carried out to execute finite element calculation on each input sample and obtain the output samples. The original data including the input and output samples are selected as the fitting samples, which are imported into Matlab. Neural network technology is used to approximate the response function with respect to the design parameters and the surrogate model is established. According to the stress-strength interface reliability model, the limit state function could be built based on the surrogate model. The reliability and reliability sensitivity analysis could be performed with the proposed copula-based high-order moment standardization technique. The specific analysis frame is shown in Figure 1.

Figure 1:

The analysis frame of the implicit limit state function problem.

5. Robust Optimization Based on Reliability Sensitivity

For the robust optimization problem of mechanical components with multiple failure modes, the reliability-based robust optimization model could be established based on the reliability optimization model, which is shown as follows:

$\begin{array}{cc}\min & f\left(\mathbf{Y}\right)={\displaystyle \sum }_{i=\mathrm{1}}^k{w}_i{f}_i\left(\mathbf{Y}\right)\\ \text{s}.\text{t}.& R_s\ge R_{\mathrm{0}}\\ & {\mathbf{Y}}_{\text{LB}}\le \mathbf{Y}\le {\mathbf{Y}}_{\text{UB}}\\ & h_i\left(\mathbf{X}\right)=\mathrm{0}\\ & q_i\left(\mathbf{X}\right)\ge \mathrm{0},\end{array}$

(38)

where f _i (·) are subobjective functions, ω _i are weighting coefficients of the subobjective functions and 0<ω _k <1, the value of which is determined by the importance degree of the subobjective function, $\mathbf{Y}={[Y_{\mathrm{1}},Y_{\mathrm{2}},\dots ,Y_n]}^{\text{T}}$ is the design variable vector, R0 is the design reliability that the mechanical component must satisfy, and Y_LB and Y_UB are the lower and upper bound, respectively. h _i (·) and q _j (·) are the equality and inequality constrains, respectively.

Two target functions are selected during the optimization process; one is the structural performance indicator, such as the mass or volume of the structure, and the other one is the reliability sensitivity which is transformed into the following form:

$f\left(\mathbf{Y}\right)={\parallel \nabla P\left(\mathbf{Y}\right)\parallel }_F=\sqrt{{\displaystyle \sum }_{i=\mathrm{1}}^k\frac{\partial R}{\partial Y_i}},$

(39)

where ∇P(Y) denotes the first-order partial derivative vector and ∥·∥F is the Frobenius norm.

Therefore, two weighting coefficients are taken into account and the corresponding expression is as follows:

$\begin{array}{c}{w}_{\mathrm{1}}=\frac{f_{\mathrm{2}}\left({\mathbf{Y}}^{*\mathrm{1}}\right)-f_{\mathrm{2}}\left({\mathbf{Y}}^{*\mathrm{2}}\right)}{\left[f_{\mathrm{1}}\left({\mathbf{Y}}^{*\mathrm{2}}\right)-f_{\mathrm{1}}\left({\mathbf{Y}}^{*\mathrm{1}}\right)\right]+\left[f_{\mathrm{2}}\left({\mathbf{Y}}^{*\mathrm{1}}\right)-f_{\mathrm{2}}\left({\mathbf{Y}}^{*\mathrm{2}}\right)\right]},\\ w_{\mathrm{2}}=\mathrm{1}-w_{\mathrm{1}}.\end{array}$

(40)

6. Numerical Examples

6.1. Reliability-Based Robust Design of an Arbor

Strength failure and stiffness failure are two common failure modes of mechanical components and are generally correlated. An arbor is taken as a study object which is simplified as a cantilever beam with strength and stiffness failure. Then the reliability, reliability sensitivity, and robust optimal design are carried out. The probabilistic properties of the random variables are listed in Table 1. F is the radial cutting force, l and d are the length and the diameter of the arbor, and E is the material elastic modulus.

Table 1:

The statistical properties of random variables of the arbor.

Random variables	Mean	Standard deviation	The third moment	The fourth moment	Distribution
F (N)	200.0	1.5	1.094	17.790	Lognormal
l (mm)	80.0	0.40	0	0.0768	Normal
d (mm)	12.0	0.60	0	0.3888	Normal
E (MPa)	2.13 × 105	5.46 × 103	5.2738 × 1010	3.123 × 1015	Lognormal

The strength failure mode and stiffness failure mode could be expressed as follows:

$\begin{array}{c}{g}_{\mathrm{1}}=\left[\sigma \right]-\frac{\mathrm{32}Fl}{\pi d^{\mathrm{3}}},\\ g_{\mathrm{2}}=\left[\omega \right]-\frac{\mathrm{64}F{l}^{\mathrm{3}}}{\mathrm{3}\pi E{d}^{\mathrm{4}}},\end{array}$

(41)

where the allowable deflection $[\omega ]=0.23$ mm and allowable stress $[\sigma ]=125$ MPa.

6.1.1. Reliability Analysis considering Correlated Failure Modes

As mentioned above, the Gumbel copula function is selected to model the correlation between failure modes. The unknown parameter in the Gumbel copula is estimated according to the data of the random variables u and v. the Kendall rank relational coefficient here is τ = 0.9081, which indicates that the two failure modes are highly correlated.

According to (14), the value of the joint probability distribution $C(P_{\mathrm{1}},P_{\mathrm{2}})=\mathrm{0.0096}$ then the failure probability and the reliability of the arbor could be obtained according to (35):

$\begin{array}{c}P=P_{\mathrm{1}}+P_{\mathrm{2}}-C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)\\ =\mathrm{0.0151}+\mathrm{0.0113}-\mathrm{0.0096}=\mathrm{0.0168},\\ R=1-P=0.9832,\\ R_{\text{MCS}}=\mathrm{0.9664},\end{array}$

(42)

where R_MCS is the reliability computed by Monte Carlo simulation.

The relative error is

$\epsilon =\frac{\left|R-R_{\text{MCS}}\right|}{R_{\text{MCS}}}=\mathrm{1.7}\%.$

(43)

6.1.2. The Reliability Sensitivity Analysis

According to (36), the reliability sensitivity matrix $\text{D}P/\text{D}\stackrel{-}{\mathbf{X}}$ and the dimensionless reliability sensitivity matrix $\text{D}{P}_w/\text{D}\stackrel{-}{\mathbf{X}}$ could be obtained as follows:

$\begin{array}{c}\frac{\text{D}P}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P}{\partial F}\\ \frac{\partial P}{\partial l}\\ \frac{\partial P}{\partial d}\\ \frac{\partial P}{\partial E}\end{array}\right]=\left[\begin{array}{c}+\mathrm{0.1031}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}\\ +\mathrm{0.2734}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}\\ -\mathrm{0.6717}\\ -\mathrm{0.2920}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\end{array}\right],\\ \frac{\text{D}{P}_w}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P}{\partial F}\\ \frac{\partial P}{\partial l}\\ \frac{\partial P}{\partial d}\\ \frac{\partial P}{\partial E}\end{array}\right]\times \frac{{\mathit{\sigma }}_{\mathbf{X}}}{P}=\left[\begin{array}{c}+\mathrm{0.0157}\\ +\mathrm{0.0111}\\ -\mathrm{0.4099}\\ -\mathrm{0.0016}\end{array}\right].\end{array}$

(44)

In this example, it is displayed that the failure probability of the arbor increases as the random variables F and l increase but descends as the random variables d and E rise. In other words, the results show that the reliability has a high dependency on d, moderate dependency on F and l, and slight dependency on E. Thus, it is necessary to reduce the reliability sensitivity of d when reliability optimization is conducted in order to reduce dependence on the random variable d.

Similarly, the reliability sensitivity matrix $\text{D}P/\text{DVar}(\mathbf{X})$ and the dimensionless matrix $\text{D}{P}_w/\text{DVar}(\mathbf{X})$ with respect to the variance of the random variable can be obtained as follows:

$\begin{array}{c}\frac{\text{D}P}{\text{DVar}\left(\mathbf{X}\right)}=\left[\begin{array}{c}\frac{\partial P}{\partial \mathrm{Var}\left(F\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(l\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(d\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(E\right)}\end{array}\right]=\left[\begin{array}{c}+\mathrm{4.6399}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{8.4227}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{1.3748}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}\\ +\mathrm{9.7380}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{12}}\end{array}\right],\\ \frac{\text{D}{P}_w}{\text{DVar}\left(\mathbf{X}\right)}=\left[\begin{array}{c}\frac{\partial P}{\partial \mathrm{Var}\left(F\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(l\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(d\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(E\right)}\end{array}\right]\times \frac{{\mathbf{\eta }}_{\mathbf{X}}}{P}=\left[\begin{array}{c}+\mathrm{1.0618}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{1.3707}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{5.0338}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\\ +\mathrm{2.9527}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\end{array}\right].\end{array}$

(45)

According to the dimensionless matrix $\text{D}{P}_w/\text{DVar}(\mathbf{X})$ , it is shown that the failure probability of the arbor increases as the variance increases. The variance of d has the greatest impact on the failure probability; then as F and l, random variable E has less impact on the failure probability.

6.1.3. Reliability-Based Robust Design

In this example, the design variable vector is defined as $\mathbf{Y}=[Y_{\mathrm{1}},Y_{\mathrm{2}}]={[l,d]}^{\text{T}}$ and the initial value vector is given as ${\mathbf{Y}}_{\mathrm{0}}={[l,d]}^{\text{T}}={[\mathrm{80,12}]}^{\text{T}}$ . According to the reliability-based robust optimal model, the optimal model of the arbor could be built as follows:

$\begin{array}{c}\min \hspace{1em}f\left(\mathbf{Y}\right)=w_{\mathrm{1}}{f}_{\mathrm{1}}\left(\mathbf{Y}\right)+w_{\mathrm{2}}{f}_{\mathrm{2}}\left(\mathbf{Y}\right)\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}=w_{\mathrm{1}}\pi \frac{l{d}^{\mathrm{2}}}{\mathrm{4}}+w_{\mathrm{2}}\sqrt{{{\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{2}}\left(\frac{\partial P}{\partial {\mathbf{Y}}_i}\right)}^{\mathrm{2}}}\\ \text{s}.\text{t}.\hspace{1em}R-R_{\mathrm{0}}>\mathrm{0.999970}\le l\le \mathrm{9010}\le d\le \mathrm{14}.\end{array}$

(46)

After optimization, the design variable vector, reliability, and reliability sensitivity become

$\begin{array}{c}{\mathbf{Y}}^{*}={\left[l,d\right]}^{\text{T}}={\left[\mathrm{70.0876,12.2548}\right]}^{\text{T}},\\ R^{*}=\mathrm{0.999984},\\ \frac{\partial P^{*}}{\partial \stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P}{\partial l}\\ \frac{\partial P}{\partial d}\end{array}\right]=\left[\begin{array}{c}+\mathrm{1.0835}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ -\mathrm{2.8397}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\end{array}\right].\end{array}$

(47)

The optimization results show that the order of magnitude of reliability sensitivity reduces a lot while meeting reliability restraint, which means that the reliability of the arbor is less dependent on the randomness of the random variables.

6.2. Reliability-Based Robust Design of a Drive Shaft

A machine drive shaft under the effect of torque is considered here. The maximum torque on the dangerous section is T, the diameter and the length of the drive shaft are d and l, the revolution is n(1500 + 180 r/min⁡), and the weight is 60 kg. The weight and the effective diameter of the circular disk are 80 kg and 400 mm. There are mainly three failure modes: the static strength failure, the fatigue strength failure, and the torsional stiffness failure. The performance functions are as follows.

The static strength failure is

$g_{\mathrm{1}}={\tau }_s-\frac{\mathrm{1.3}T}{W_{\tau }}.$

(48)

The fatigue strength failure is

$g_{\mathrm{2}}={\tau }_{\mathrm{0}K}-\frac{T}{W_{\tau }}.$

(49)

The torsional stiffness failure is

$g_{\mathrm{3}}=\left[\phi \right]-\frac{Tl}{GI},$

(50)

where τ _s is the torsion yield limit of the drive shaft and τ _s = 0.5σ _s , σ _s is the yield limit the and σ _s = 280 MPa, T is the torque, W_τ is the shear section, and $W_{\tau }=\pi d^{\mathrm{3}}/\mathrm{16}$ . τ_0K is the pulsating cycle fatigue limit, τ₋₁ is the torsion fatigue limit and ${\tau }_{\mathrm{0}K}=\mathrm{1.5}\hspace{0.17em}{\tau }_{-\mathrm{1}}/\mathrm{1.904}=\mathrm{0.1229}({\sigma }_s+{\sigma }_b)$ , σ _b is the material strength limit and σ _b = 610 MPa, $[\phi ]$ is the torsional stiffness allowance and $[\phi ]=\mathrm{0.0105}$ , G is shear modulus of elasticity, and I is the moment of inertia. Here σ _s , σ _b , d, T, and l are considered as the random variables and the corresponding statistic properties are listed in Table 2.

Table 2:

The statistical properties of random variables of the drive shaft.

Random variables	Mean	Standard deviation	Distribution
d (mm)	150	0.75	Normal
l (mm)	600	3.00	Normal
T (kmm)	6.00 × 107	6.00 × 106	Normal
σ s (MPa)	320	32.00	Normal
σ b (MPa)	610	30.50	Normal

6.2.1. Reliability Analysis considering Correlated Failure Modes

Based on the known probability information, the reliability index, the reliability, and the failure probability of each failure mode could be calculated by the proposed method:

$\begin{array}{c}{\beta }_{\mathrm{1}}=\mathrm{2.1210},\hspace{1em}\hspace{1em}{R}_{\mathrm{1}}=\mathrm{0.9830},\hspace{1em}\hspace{1em}{P}_{\mathrm{1}}=\mathrm{0.0170},\\ {\beta }_{\mathrm{2}}=\mathrm{2.2313},\hspace{1em}\hspace{1em}{R}_{\mathrm{2}}=\mathrm{0.9872},\hspace{1em}\hspace{1em}{P}_{\mathrm{2}}=\mathrm{0.0128},\\ {\beta }_{\mathrm{3}}=\mathrm{1.6646},\hspace{1em}\hspace{1em}{R}_{\mathrm{3}}=\mathrm{0.9520},\hspace{1em}\hspace{1em}{P}_{\mathrm{3}}=\mathrm{0.0480}.\end{array}$

(51)

According to the Gumbel function, the $C(P_i,P_h)\hspace{0.17em}\hspace{0.17em}(\mathrm{1}\le i,\hspace{0.17em}h\le \mathrm{3})$ and $C(P_i,P_h,P_t)\hspace{0.17em}\hspace{0.17em}(\mathrm{1}\le i<h<t\le \mathrm{3})$ could be obtained on the basis of the above results:

$\begin{array}{c}C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)=\mathrm{0.0088},\\ C\left(P_{\mathrm{1}},P_{\mathrm{3}}\right)=\mathrm{0.0126},\\ C\left(P_{\mathrm{2}},P_{\mathrm{3}}\right)=\mathrm{0.0121},\\ C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}\right)=C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),P_{\mathrm{3}}\right)=\mathrm{0.0085}.\end{array}$

(52)

The spearman rank correlation coefficients between two failure modes could then be obtained as follows, which show that every two failure modes are highly correlated:

$\begin{array}{c}{\tau }_{\mathrm{12}}=\mathrm{0.9912},\\ {\tau }_{\mathrm{13}}=\mathrm{0.9874},\\ {\tau }_{\mathrm{23}}=\mathrm{0.9879}.\end{array}$

(53)

According to the copula theory, the failure probability P and the reliability R can be obtained while considering three correlated failure modes:

$\begin{array}{c}P={\displaystyle \sum }_{i=\mathrm{1}}^N{P}_i-{\displaystyle \sum }_{\mathrm{1}\le i<h\le \mathrm{3}}C\left(P_i,P_h\right)+{\displaystyle \sum }_{\mathrm{1}\le i<h<t\le \mathrm{3}}C\left(P_i,P_h,P_t\right)\\ =P_{\mathrm{1}}+P_{\mathrm{2}}+P_{\mathrm{3}}-C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)-C\left(P_{\mathrm{1}},P_{\mathrm{3}}\right)\\ \hspace{1em}-C\left(P_{\mathrm{2}},P_{\mathrm{3}}\right)+C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}\right)=\mathrm{0.0527},\\ R=\mathrm{1}-P=\mathrm{0.9473}.\end{array}$

(54)

The reliability calculated by Monte Carlo simulation with 10⁴ samples is

$R_{\text{MCS}}=\mathrm{0.9460}.$

(55)

The relative error is

$\epsilon =\frac{\left|R-R_{\text{MCS}}\right|}{R_{\text{MCS}}}=\mathrm{0.14}\%.$

(56)

6.2.2. The Reliability Sensitivity Analysis

According to (36), the reliability sensitivity analysis could then be performed. The reliability sensitivity matrix $\text{D}P/\text{D}\stackrel{-}{\mathbf{X}}$ and the dimensionless matrix $\text{D}{P}_w/\text{D}\stackrel{-}{\mathbf{X}}$ with respect to the mean of the random variable can be obtained as follows:

$\begin{array}{c}\frac{\text{D}P}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P}{\partial {\sigma }_s}\\ \frac{\partial P}{\partial {\sigma }_b}\\ \frac{\partial P}{\partial T}\\ \frac{\partial P}{\partial d}\\ \frac{\partial P}{\partial l}\end{array}\right]=\left[\begin{array}{c}-\mathrm{4.2148}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ -\mathrm{3.0485}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{1.6785}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{8}}\\ -\mathrm{2.2900}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ +\mathrm{1.4877}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\end{array}\right],\\ \frac{\text{D}{P}_w}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P}{\partial {\sigma }_s}\\ \frac{\partial P}{\partial {\sigma }_b}\\ \frac{\partial P}{\partial T}\\ \frac{\partial P}{\partial d}\\ \frac{\partial P}{\partial l}\end{array}\right]\times \frac{{\mathit{\sigma }}_{\mathbf{X}}}{P}=\left[\begin{array}{c}-\mathrm{0.2559}\\ -\mathrm{0.0176}\\ +\mathrm{1.9110}\\ -\mathrm{0.0326}\\ +\mathrm{0.0847}\end{array}\right].\end{array}$

(57)

According to the dimensionless matrix $\partial P_w/\partial {\stackrel{-}{\mathbf{X}}}^{\text{T}}$ , it is shown that the reliability of the drive shaft increases as random variables σ _s , σ _b , and d increase. On the contrary, while random variables T and l increase, the reliability of the drive shaft decreases. In other words, the mean of T has the greatest impact on the structural reliability; then as σ _s and σ _b , random variables d and l have less impact on the failure probability. Based on the analysis, it is necessary to control the maximum torque T in order to guarantee the target reliability of the drive shaft under certain load cases.

Similarly, the reliability sensitivity matrix $\text{D}P/\text{DVar}(\mathbf{X})$ and the dimensionless matrix $\text{D}{P}_w/\text{DVar}(\mathbf{X})$ with respect to the variance of the random variable can be obtained as follows:

$\begin{array}{c}\frac{\text{D}P}{\text{DVar}\left(\mathbf{X}\right)}=\left[\begin{array}{c}\frac{\partial P}{\partial \mathrm{Var}\left({\sigma }_s\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left({\sigma }_b\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(T\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(d\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(l\right)}\end{array}\right]=\left[\begin{array}{c}+\mathrm{1.0789}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{3.9262}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{7}}\\ +\mathrm{2.2803}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{15}}\\ +\mathrm{3.1569}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{2.0610}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\end{array}\right],\\ \frac{\text{D}{P}_w}{\text{DVar}\left(\mathbf{X}\right)}=\left[\begin{array}{c}\frac{\partial P}{\partial \mathrm{Var}\left({\sigma }_s\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left({\sigma }_b\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(T\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(d\right)}\\ \frac{\partial P}{\partial \mathrm{Var}\left(l\right)}\end{array}\right]\times \frac{{\mathbf{\eta }}_{\mathbf{X}}}{P}=\left[\begin{array}{c}+\mathrm{0.2096}\\ +\mathrm{0.0069}\\ +\mathrm{1.5577}\\ +\mathrm{0.0034}\\ +\mathrm{0.0035}\end{array}\right].\end{array}$

(58)

It can be seen from the reliability sensitivity matrix that the variance of random variables has negative impact on the reliability of the drive shaft, which means that the increase of the variance would cause the component tend to fail. Among the design variables, the torque T has the greatest impact on the structural reliability; then as σ _s , random variables d, l, and σ _b have less impact on the failure probability.

6.2.3. The Reliability-Based Robust Design

Based on the results of the reliability and reliability sensitivity, the reliability-based robust optimal design could be conducted. The geometric variables of the drive shaft are selected as the design variables, being $\mathbf{Y}={[d,l]}^{\text{T}}$ . The model of the optimization could then be built according to the reliability robust optimization theory:

$\begin{array}{c}\min \hspace{1em}\hspace{0.17em}f=w_{\mathrm{1}}{f}_{\mathrm{1}}+w_{\mathrm{2}}{f}_{\mathrm{2}}=w_{\mathrm{1}}\pi {\left(\frac{d}{\mathrm{2}}\right)}^{\mathrm{2}}l+w_{\mathrm{2}}{\parallel \nabla P\left(\mathbf{Y}\right)\parallel }_F\\ \text{s}.\text{t}.\hspace{1em}\hspace{1em}R-R_{\mathrm{0}}\ge \mathrm{0130}\le d\le \mathrm{170500}\le l\le \mathrm{610}.\end{array}$

(59)

The subobjective function $f_{\mathrm{1}}=\pi {(d/\mathrm{2})}^{\mathrm{2}}l$ represents the volume of the drive shaft; another subobjective function f₂ = ∥∇P(Y)∥2 is the F-norm of the reliability sensitivity of the design variables. ∇P(Y) is the vector of the reliability sensitivity with respect to the mean. R0 is the reliability target and R₀ = 0.99. The initial value of the design variable vector is ${\mathbf{Y}}_{\mathrm{0}}={[\mathrm{150,530}]}^{\text{T}}$ .

After optimization, the optimal solution of the design variable vector is ${\mathbf{Y}}^{*}={[\mathrm{153.22,500.99}]}^{\text{T}}$ . The reliability and the reliability sensitivity of the drive shaft after optimization could then be obtained as follows:

$\begin{array}{c}{R}^{*}=\mathrm{0.99395},\\ \frac{\text{D}{P}^{*}}{\text{D}\stackrel{-}{\mathbf{Y}}}=\left[\begin{array}{c}\frac{\partial P^{*}}{\partial d}\\ \frac{\partial P^{*}}{\partial l}\end{array}\right]=\left[\begin{array}{c}-\mathrm{7.3204}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{2.7220}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\end{array}\right].\end{array}$

(60)

The results obtained show that the reliability of the structure after optimization meets the target reliability requirement, and the reliability sensitivity is also decreased. The reliability of the optimized drive shaft is less dependent on the randomness of the random variables.

6.3. The Implicit Limit State Function Problem: Reliability-Based Robust Design of a Bracket

As shown in Figure 2 is a bracket which is parameterized with some critical structure parameters. In this example, three failure modes of the bracket, namely, the strength failure under vertical load, the displacement excursion, and the strength failure under combined loads, are defined. The proposed optimal method for implicit limit state function problem as mentioned above is applied to put forward optimization scheme for the bracket.

Figure 2:

The 3-dimensional structure of the bracket.

The deterministic analysis of the bracket by FEM is performed and one of the results is illustrated in Figure 3, from which the location and the value of the maximum stress are clearly shown.

Figure 3:

The stress nephogram of the bracket.

According to the modeling method of the surrogate model, the DOE is firstly carried out. The input samples are generated by Latin hypercube sampling (LHS) method from the input variables. Then, the input sample matrix (N × M) is established in iSIGHT, (N represents the repeat number of the DOE, and M represents the number of the design variables). During each experiment, the model is rebuilt through dimension-driven method in Pro/E, and the response is analyzed in MSC. Nastran.

The data of the input-output samples could be obtained by performing DOE study. The number of the samples depends on the modeling method of the surrogate model and a test sample which could test the generalization ability of the established surrogate model is always recommended. Here the ANN method is adopted to establish the surrogate model of the response. The combined load condition of the bracket is taken as the example, and the relative error of the fitting function is shown in Figure 4.

Figure 4:

The relative error of the fitting function.

Based on the approximated function of the response, the limit state function for reliability analysis could be further built. Let the response function be S = g(X) and the generalized allowable stress Q; then the limit state function could be expressed as

$Z=Q-S=Q-g\left(\mathbf{X}\right),$

(61)

where $\mathbf{X}={[f_{\mathrm{1}},f_{\mathrm{2}},f_{\mathrm{3}},L_{\mathrm{1}},L_{\mathrm{2}},L_{\mathrm{3}},L_{\mathrm{4}}]}^{\text{T}}$ represents the random variable vector, the probabilistic properties of which are shown in Table 3.

Table 3:

The probabilistic properties of random variables.

Random variables	Mean	Standard variance	The third moment	The fourth moment	Distribution
f1 (N)	1470	7.35	1.287 × 102	1.026 × 104	Lognormal
f2 (N)	1100	22.0	3.450 × 103	8.232 × 105	Lognormal
f3 (N)	1000	20.0	2.592 × 103	5.622 × 105	Lognormal
L1 (mm)	200	1.0	0	3.0	Normal
L2 (mm)	700	3.5	0	450.188	Normal
L3 (mm)	200	1.0	0	3.0	Normal
L4 (mm)	120	0.6	0	0.389	Normal

6.3.1. Reliability and Reliability Sensitivity Analysis under Each Failure Mode

The stress analysis of the bracket subjected to vertical load is firstly carried out and the failure mode one is defined as the maximum stress exceeding the stress limitation; then the failure probability and reliability could be calculated as follows:

$P_{\mathrm{1}}=\mathrm{0.0012},\hspace{1em}\hspace{1em}{R}_{\mathrm{1}}=\mathrm{0.9988}.$

(62)

The design variable vector is ${\mathbf{X}}_{\mathrm{1}}={[f_{\mathrm{1}},L_{\mathrm{1}},L_{\mathrm{2}},L_{\mathrm{3}},L_{\mathrm{4}}]}^{\text{T}}$ under this failure mode; according to the reliability sensitivity analysis method mentioned above, the reliability sensitivity $\text{D}{P}_{\mathrm{1}}/\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}$ and the dimensionless transition $\text{D}{P}_{\mathrm{1}w}/\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}$ could be obtained as

$\begin{array}{c}\frac{\text{D}{P}_{\mathrm{1}}}{\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}}=\left[\begin{array}{c}\frac{\partial P_{\mathrm{3}}}{\partial f_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{4}}}\end{array}\right]=\left[\begin{array}{c}+\mathrm{1.1943}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\\ +\mathrm{3.6736}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{1.7021}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\\ -\mathrm{1.8326}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ -\mathrm{2.7870}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\end{array}\right],\\ \frac{\text{D}{P}_{\mathrm{1}w}}{\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}}=\left[\begin{array}{c}\frac{\partial P_{\mathrm{3}}}{\partial f_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{4}}}\end{array}\right]\times \frac{{\mathit{\sigma }}_{{\mathbf{X}}_{\mathrm{1}}}}{P_{\mathrm{1}}}=\left[\begin{array}{c}+\mathrm{7.2865}\times \mathrm{1}{\mathrm{0}}^{\mathrm{1}}\\ +\mathrm{3.0494}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\\ +\mathrm{4.9451}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ -\mathrm{1.5213}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\\ +\mathrm{1.3881}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\end{array}\right].\end{array}$

(63)

The failure mode two is defined as the displacement of the bracket exceeding the displacement limitation; then the failure probability and reliability of the bracket could be calculated as follows:

$P_{\mathrm{2}}=\mathrm{0.0212},\hspace{1em}\hspace{1em}{R}_{\mathrm{2}}=\mathrm{0.9788}.$

(64)

The reliability sensitivity $\text{D}{P}_{\mathrm{2}}/\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}$ and the dimensionless transition $\text{D}{P}_{\mathrm{2}w}/\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}$ could be obtained as

$\begin{array}{c}\frac{\text{D}{P}_{\mathrm{2}}}{\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}}=\left[\begin{array}{c}\frac{\partial P_{\mathrm{2}}}{\partial f_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{4}}}\end{array}\right]=\left[\begin{array}{c}+\mathrm{1.5746}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}\\ -\mathrm{5.6819}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{1.8122}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{3.8587}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{1.3139}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\end{array}\right],\\ \frac{\text{D}{P}_{\mathrm{2}w}}{\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{1}}}=\left[\begin{array}{c}\frac{\partial P_{\mathrm{2}}}{\partial f_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{2}}}{\partial L_{\mathrm{4}}}\end{array}\right]\times \frac{{\mathit{\sigma }}_{{\mathbf{X}}_{\mathrm{1}}}}{P_{\mathrm{2}}}=\left[\begin{array}{c}+\mathrm{5.4521}\times \mathrm{1}{\mathrm{0}}^{\mathrm{1}}\\ -\mathrm{2.6767}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ +\mathrm{2.9879}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ +\mathrm{1.8178}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ +\mathrm{3.7139}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\end{array}\right].\end{array}$

(65)

Similarly, failure mode three is defined as the maximum stress exceeding the stress limitation, when the bracket is subjected to the vertical load f2 and also the lateral load f3. The failure probability and reliability under this failure mode could be obtained as follows:

$P_{\mathrm{3}}=\mathrm{0.0011},\hspace{1em}\hspace{1em}{R}_{\mathrm{3}}=\mathrm{0.9989}.$

(66)

In failure mode three, the design variable vector is ${\mathbf{X}}_{\mathrm{2}}={[f_{\mathrm{2}},f_{\mathrm{3}},L_{\mathrm{1}},L_{\mathrm{2}},L_{\mathrm{3}},L_{\mathrm{4}}]}^{\text{T}}$ and the corresponding reliability sensitivity matrix could be obtained as

$\begin{array}{c}\frac{\text{D}{P}_{\mathrm{3}}}{\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{2}}}=\left[\begin{array}{c}\frac{\partial P_{\mathrm{3}}}{\partial f_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial f_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{4}}}\end{array}\right]=\left[\begin{array}{c}+\mathrm{1.7528}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{1.4388}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{2.2146}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\\ -\mathrm{1.1192}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\\ +\mathrm{1.2341}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ -\mathrm{1.2209}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\end{array}\right],\\ \frac{\text{D}{P}_{\mathrm{3}w}}{\text{D}{\stackrel{-}{\mathbf{X}}}_{\mathrm{2}}}=\left[\begin{array}{c}\frac{\partial P_{\mathrm{3}}}{\partial f_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial f_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_{\mathrm{3}}}{\partial L_{\mathrm{4}}}\end{array}\right]\times \frac{{\mathit{\sigma }}_{{\mathbf{X}}_{\mathrm{2}}}}{P_{\mathrm{3}}}=\left[\begin{array}{c}+\mathrm{3.5552}\\ +\mathrm{2.6529}\\ +\mathrm{2.0417}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ -\mathrm{3.6114}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\\ +\mathrm{1.1378}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\\ -\mathrm{6.7535}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{3}}\end{array}\right].\end{array}$

(67)

6.3.2. Reliability and Reliability Sensitivity Analysis of the Bracket with Correlated Failure Modes

The three failure modes are defined as above, and the failure probability and reliability under each failure mode are calculated, respectively. However, the three potential failure modes are not independent; one failure of the bracket may lead to another failure to some extent. According to the proposed copula-based method, the Gumbel copula function is used to describe the correlation between each failure mode. The value of the joint distribution function of failure modes could be obtained as follows:

$\begin{array}{c}C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right)=\mathrm{0.0010},\\ C\left(P_{\mathrm{1}},P_{\mathrm{3}}\right)=\mathrm{0.0010},\\ C\left(P_{\mathrm{2}},P_{\mathrm{3}}\right)=\mathrm{9.2891}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}},\\ C\left(P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}\right)=C\left(C\left(P_{\mathrm{1}},P_{\mathrm{2}}\right),P_{\mathrm{3}}\right)=\mathrm{2.8672}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}.\end{array}$

(68)

Based on the obtained results, the Spearman rank correlation coefficients of any two failure modes are

$\begin{array}{c}{\tau }_{\mathrm{12}}=\mathrm{0.9820},\\ {\tau }_{\mathrm{13}}=\mathrm{0.7490},\\ {\tau }_{\mathrm{23}}=\mathrm{0.7495}.\end{array}$

(69)

As can be seen from the correlation coefficient, the failure modes one and two are highly correlated; the failure modes one and three and the failure modes two and three also have a strong correlation. While considering the impact of the correlation, the failure probability, reliability, and the reliability sensitivity of the bracket could be calculated:

$\begin{array}{c}{P}_S=\mathrm{0.0208},\\ R_S=\mathrm{0.9792},\\ \frac{\text{D}{P}_S}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P_S}{\partial f_{\mathrm{1}}}\\ \frac{\partial P_S}{\partial f_{\mathrm{2}}}\\ \frac{\partial P_S}{\partial f_{\mathrm{3}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{4}}}\end{array}\right]=\left[\begin{array}{c}+\mathrm{3.7523}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{2}}\\ +\mathrm{1.7257}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{1.4165}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ -\mathrm{2.8061}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\\ +\mathrm{3.2849}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{6}}\\ +\mathrm{1.6916}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\\ +\mathrm{1.1407}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{5}}\end{array}\right].\end{array}$

(70)

The dimensionless transition is

$\frac{\text{D}{P}_{Sw}}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P_S}{\partial f_{\mathrm{1}}}\\ \frac{\partial P_S}{\partial f_{\mathrm{2}}}\\ \frac{\partial P_S}{\partial f_{\mathrm{3}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_S}{\partial L_{\mathrm{4}}}\end{array}\right]\times \frac{{\mathit{\sigma }}_{\mathbf{X}}}{P_S}=\left[\begin{array}{c}+\mathrm{1.3257}\times \mathrm{10}\\ +\mathrm{1.8248}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}\\ +\mathrm{1.3617}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{1}}\\ -\mathrm{1.3488}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{5.5264}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{8.1311}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\\ +\mathrm{3.2897}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{4}}\end{array}\right].$

(71)

According to the results, when considering the correlation of the failure modes, it shows that the failure probability of the bracket increases as the load (f₁,f₂,f₃) and geometry (L₂,L₃,L₄) increase but descends as the geometry L1 rises. Besides, in terms of the degree of influences, the reliability is very sensitive to the load f1, moderately sensitive to the loads f2 and f3, and slightly sensitive to the geometries.

6.3.3. Reliability-Based Robust Design of the Bracket with Correlated Failure Modes

Based on the reliability and reliability sensitivity results, the reliability-based optimization could be further performed, the model of the reliability sensitivity-based robust design can be built as

$\begin{array}{c}\min \hspace{1em}f\left(\mathbf{Y}\right)=w_{\mathrm{1}}{f}_{\mathrm{1}}\left(\mathbf{Y}\right)+w_{\mathrm{2}}{f}_{\mathrm{2}}\left(\mathbf{Y}\right)\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\\ =w_{\mathrm{1}}V+w_{\mathrm{2}}\sqrt{{{\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{4}}\left(\frac{\partial P_S}{\partial {\mathbf{Y}}_i}\right)}^{\mathrm{2}}}\\ \text{s}.\text{t}.\hspace{1em}{R}_S-R_{\mathrm{0}}>\mathrm{0.9999190}\le L_{\mathrm{1}}\le \mathrm{210690}\le L_{\mathrm{2}}\le \mathrm{710190}\le L_{\mathrm{3}}\le \mathrm{210110}\le L_{\mathrm{4}}\le \mathrm{130}.\end{array}$

(72)

The design variable vector is $\mathit{Y}={[L_1,L_2,L_3,L_4]}^{\text{T}}$ with the mean value ${\mathit{Y}}_{\mathit{0}}={[\mathrm{200,700,200,120}]}^{\text{T}}$ and the target reliability is set as R₀ = 0.9999. The subobjective function $f_{\mathrm{1}}(\mathbf{Y})=V$ is the total volume of the bracket, and the subobjective function $f_{\mathrm{2}}(\mathbf{Y})=\sqrt{{{\sum }_{i=\mathrm{1}}^{\mathrm{4}}(\partial P_S/\partial {\mathbf{Y}}_i)}^{\mathrm{2}}}$ is the Frobenius norm of the reliability sensitivity with respect to the design variables, which means the overall effect of the random fluctuation of the design variables. The size constraints of the optimization problem are also defined in the model.

After optimization, the design variable vector becomes

$Y^{*}={\left[\mathrm{209.9752,690.0137,190.3129,110.6802}\right]}^{\text{T}}.$

(73)

The corresponding failure probability, reliability, and the reliability sensitivity become

$\begin{array}{c}{P}_S^{*}=\mathrm{0.00004718},\hspace{1em}\hspace{1em}{R}_S^{*}=\mathrm{0.99995282},\\ \frac{\text{D}{P}_S^{*}}{\text{D}\stackrel{-}{\mathbf{X}}}=\left[\begin{array}{c}\frac{\partial P_S^{*}}{\partial L_{\mathrm{1}}}\\ \frac{\partial P_S^{*}}{\partial L_{\mathrm{2}}}\\ \frac{\partial P_S^{*}}{\partial L_{\mathrm{3}}}\\ \frac{\partial P_S^{*}}{\partial L_{\mathrm{4}}}\end{array}\right]=\left[\begin{array}{c}+\mathrm{4.7213}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{8}}\\ +\mathrm{5.3061}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{9}}\\ +\mathrm{1.0956}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{8}}\\ +\mathrm{2.7765}\times \mathrm{1}{\mathrm{0}}^{-\mathrm{8}}\end{array}\right].\end{array}$

(74)

The optimization results show that the order of magnitude of reliability sensitivity reduces a lot while meeting reliability target and structure geometry restraint, which means that the reliability of the bracket is less dependent on the design variables, the changes of the design variables lead to less fluctuation of the reliability, and the structure reliability has robust features.

7. Conclusions

The problem of the reliability-based robust optimization of mechanical components with multiple failure modes is addressed, considering the positive correlation between failure modes. In this work, the reliability analysis is performed with the moment-based method, which is suitable for arbitrary distributed variables. The copula function theory is introduced and adopted to model the nonlinear positive correlation between failure modes. The system reliability model of the multiple failure modes is built with the theory of probability based on the copula-based joint distribution function. After solving a matrix derivation problem, the reliability sensitivity with respect to the mean and the variance of the random variables is obtained for the case of independent and correlated failure modes. Then the robust optimization of the system reliability problem is performed as a biobjective optimization problem by minimizing the reliability sensitivity and parameter constraints. For implicit limit state function problem, the design of experiment and the surrogate model are adopted to establish the response functions of failure modes, and the robust optimization is finally conducted with the proposed numerical models. The results of the numerical examples serve as testimony to the effectiveness of the proposed method and should be useful when considering the reliability problem of mechanical component with correlated failure modes.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.