Reliability analysis for high-lift device based on Copula function and Bayesian updating

Abstract

In order to improve the reliability analysis accuracy of the aircraft high-lift, an approach based on the Copula function theory and Bayesian updating is proposed. Considering the influence of the random variables’ correlation in the process of updating, choosing the reasonable prior joint distribution and likelihood function is crucial. Under the condition of the incomplete probability information, the analytic expressions of the prior joint distribution and likelihood function of the correlated random variables are derived through the Copula function. Then, the posterior joint distribution is obtained by Bayesian updating. The reliability of the lifting device is calculated based on the posterior distribution. The case analysis shows that the reliability results based on the proposed approach are more accurate and more coincident with the factual situation than the reliability analysis results based on the independence assumption of random variables.

Introduction

The high-lift device is used to improve take-off weight and add lift for take-off and landing of civil airplanes. ¹ The safety of airplanes has close relationship with the reliability of the high-lift device. With the development and service of civil airplanes around the world, airworthiness of civil airplanes puts forward higher request to the reliability of the high-lift device. Therefore, analyzing the reliability of the high-lift device accurately is a vital issue in the design phase of civil airplanes.

Focusing on the above issue, the statistical approach has been applied to calculate and analyze the reliability and failure probability of the high-lift device mechanism in the design process of civil airplanes. ² Because the sample data available are often limited in practical engineering, the probability distribution functions of random variables cannot be obtained precisely. The application of statistical methods is restricted. In order to solve the problem of data deficiency, Bayesian updating approach is adopted to obtain the more accurate probability distribution estimation of random variables through the integration of the experimental data. ³ Bayesian updating has been widely used in many engineering fields, such as machinery, electronic, and civil engineering.^4–9 Due to the computation complexity, few studies consider the correlation of random variables in the updating process. In light of this situation that the correlation has to be considered, some approximate methods are adopted to deal with the problem of the correlation of random variables. Zhu and Frangopol ³ assume that the correlated random variables follow the multidimensional normal distribution when assessing the reliability of ship structures using Bayesian updating. An et al. ¹⁰ assumed that the correlated random variables follow the multidimensional normal and lognormal distributions, and use Bayesian statistics method to estimate the wear coefficient and predict the wear volume in the revolute joint. For the sake of convenience, Lin and Chen ¹¹ adopt the assumption of independence in their research. However, in the above approximate methods, there exist obvious disadvantages: the assumption of independence may lead to large errors in some cases; different types of multidimensional joint distribution assumptions can cause different Bayesian updating results. Tang and colleagues^12,13 state that the correlation between the shear strength parameters has a great influence on the stability reliability of the slope, and the reliability of the slope will be underestimated obviously if the correlation is neglected. Tang et al. ¹⁴ also state that the simulated load–displacement hyperbolic curves of the single pile dramatically deviates from the measured load–displacement hyperbolic curves without considering the correlation between geomechanical parameters. The analysis results of An et al. ¹⁰ show that using different types of joint distributions to construct the likelihood functions leads to significant difference in the posterior distribution. In comparison, the incomplete probability information, namely, the correlation coefficients and the marginal probability distributions of random variables, is more easily obtained in engineering.^15,16 However, if the joint probability distribution function of the variables is not the multidimensional normal distribution, the prior joint distribution and the likelihood function cannot be uniquely determined by the incomplete probability information. ¹⁷ Under the condition of the incomplete probability information, how to choose a reasonable prior joint distribution and likelihood function is the key problem to be solved.

Under the condition of the incomplete probability information, the Copula function provides an effective way to construct the joint distribution of random variables. The application of Copula function theory is first seen in the data analysis of the financial finance field, ¹⁸ and there have also been research achievements about the hydrology and the civil engineering reliability in recent years.^{12–14,17,19–21} However, the research results in the field of mechanism reliability are rarely published. For the correlated random variables, the Copula function is introduced to construct the prior joint distribution and the likelihood function before Bayesian updating of the aircraft lift device. Then, the posterior joint distribution of random variables is obtained by updating, and the reliability of the lifting device is calculated.

Basic theory of Bayesian updating and Copula function

Bayesian updating

Bayesian updating theorem

Based on the traditional statistical approach, a large amount of samples are required to construct the probability distributions of random variables precisely. However, the sample data are shortage in most engineering cases. The traditional statistical approach is no longer applicable. By incorporating the prior knowledge with the new information, such as the experimental data and empirical data, we can obtain the more precise probability distribution estimations. This is the main idea of the Bayesian updating approach. ⁶

Let the symbol X be a random vector, and the prior joint probability density function (PDF) of the random vector be $f(\mathit{X})$ . Then, let the symbol Y represent the experimental data about $\mathit{X}$ . With the experimental data, we can construct the likelihood function of the random vector $\mathit{X}$ . Let the symbol $L(\mathit{Y}|\mathit{X})$ be the likelihood function of $\mathit{X}$ , and $f^{*}(\mathit{X}|\mathit{Y})$ be the posterior joint PDF of $\mathit{X}$ . The principle of the updating method is as follows ⁵

$$f^{*}\left(\mathit{X}|\mathit{Y}\right)=\mathrm{kL}\left(\mathit{Y}|\mathit{X}\right)f\left(\mathit{X}\right)$$

(1)

where k is the normalizing constant and given by ⁵

$$k={\left(\underset{-\infty }{\overset{+\infty }{\int }}L\left(\mathit{Y}|\mathit{X}\right)f\left(\mathit{X}\right)\mathrm{d}\mathit{X}\right)}^{-1}$$

(2)

By the Bayesian updating approach, the prior joint PDF of the random vector X can be updated to the posterior $f^{*}(\mathit{X}|\mathit{Y})$ . Comparing the prior $f(\mathit{X})$ with the posterior $f^{*}(\mathit{X}|\mathit{Y})$ , the addition of the experimental data $\mathit{Y}$ during the updating process causes the change of the probability distributions model of $\mathit{X}$ . In other words, the knowledge about the random vector X is updated by the Bayesian updating approach. As more information about X is obtained, the updated probability distributions of X are more accurate.

Markov chain Monte Carlo method

Through the Bayesian updating approach, we can obtain the expression of the posterior joint PDF of $\mathit{X}$ . However, when we calculate the posterior joint PDF, one of the practical problems faced by us is that computing integrals of the posterior joint PDF is very difficult, or even impossible. In this case, the numerical calculation methods become the practical choice. The Markov chain Monte Carlo (MCMC) method is the most popular and effective technique among the others in recent years and has gotten extensive application.^10,22 The MCMC method contains a number of sampling algorithms, such as the Gibbs sampling algorithm and the M-H sampling algorithm.^22,23 The M-H algorithm is chosen in this article. The simulation process of the MCMC method using the M-H algorithm is as follows:

• Step 1. Initialize the vector ${\mathit{X}}^{\left(0\right)}$ .

• Step 2. For $i=0\mathrm{to}N-1$ :

Sample $u$ ∼ $\mathrm{U}(0,1)$

• Step 3.

If

$$u<min(1,\frac{f^{*}\left({\mathit{X}}^{*}\right)q\left({\mathit{X}}^{\left(i\right)}|{\mathit{X}}^{*}\right)}{f^{*}\left({\mathit{X}}^{\left(i\right)}\right)q\left({\mathit{X}}^{*}|{\mathit{X}}^{\left(i\right)}\right)})$$

Then

$${\mathit{X}}^{(i+1)}={\mathit{X}}^{*}$$

Else

$${\mathit{X}}^{(i+1)}={\mathit{X}}^{\left(i\right)}$$

where ${\mathit{X}}^{\left(0\right)}$ is the initial values of the random vector $\mathit{X}$ , N is the number of samples, $\mathrm{U}(0,1)$ is the uniform distribution in the interval of [0, 1], and $q(\mathit{X})$ is the proposal probability function or the transition function.

When there are a large enough number of samples, we can well get the posterior distribution characteristics of the random vector. According to previous studies, some approaches have been provided to diagnose and analyze the convergence of the MCMC method. ²³ For reasons of simplicity, the graphical method is adopted to analyze the convergence in this article.

Copula function

In the Bayesian updating process, ignoring the correlation of random variables will lead to imprecise posterior probability distribution and further lead to even wrong reliability analysis results. Therefore, constructing the prior joint PDF and likelihood function of the correlated random variables is crucial. Due to the limited data, we can only obtain the marginal probability distributions and correlation coefficient.^15,16 The issue is that the prior joint PDF and likelihood function cannot be determined uniquely by these data. ¹⁷ To solve the issue, the Copula function theory is applied in this article. With the Copula function, we can build the prior joint PDF and likelihood function of correlated random variables using only the marginal probability distributions and correlation coefficient. Sklar first proposed the Copula function theory in 1959. ²⁴ The Sklar’s theorem is as follows.

Assume the symbol $F(x_1,x_2,\dots ,x_n)$ represents the joint cumulative distribution function (CDF) of the correlated random variables $x_1,x_2,\dots ,x_n$ , and the marginal CDFs of them are, respectively, $F_1(x_1),F_2(x_2),\dots ,F_n(x_n)$ . Then, there is an n-dimensional Copula function $C(·)$ as follows

$$F(x_1,x_2,\dots ,x_n)=C(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_n\left(x_n\right);\theta )$$

(3)

where $\theta$ is the Copula parameter which is used to express the correlation of the correlated random variables.

If the function $C(·)$ is a Copula function, and at the same time, functions $F_1(x_1),F_2(x_2),\dots ,F_n(x_n)$ are the CDFs of random variables $x_1,x_2,\dots ,x_n$ , respectively, the function $F(x_1,x_2,\dots ,x_n)$ in equation (3) is the joint CDF of these random variables. ²⁴ Especially, the marginal CDFs of $F(x_1,x_2,\dots ,x_n)$ are $F_1(x_1),F_2(x_2),\dots ,F_n(x_n)$ . Based on the definition of Sklar’s theorem, we can get the prior joint PDF and likelihood function of the correlated random variables by simply knowing a Copula function and their marginal distributions. This is especially significant in the practical engineering application.

Furthermore, through calculating the derivative of $F(x_1,x_2,\cdots ,x_n)$ , we can obtain the joint PDF of $x_1,x_2,\dots ,x_n$ , which is given by

$$\begin{array}{c}\multicolumn{1}{c}{f(x_1,x_2,\dots ,x_n)}\\ \multicolumn{1}{c}{=c(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_n\left(x_n\right);\theta )f_1\left(x_1\right)f_2\left(x_2\right)\cdots f_n\left(x_n\right)}\end{array}$$

(4)

in which $f_1(x_1),f_2(x_2),\dots ,f_n(x_n)$ are the marginal PDFs of $x_1,x_2,\dots ,x_n$ , respectively, and $c(F_1(x_1),F_2(x_2),\dots ,F_n(x_n);\theta )$ represents the Copula density function of $x_1,x_2,\dots ,x_n$ . Through calculating the derivative of the Copula function $C(·)$ , we can acquire the Copula density function

$$\begin{array}{c}\multicolumn{1}{c}{c(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_n\left(x_n\right);\theta )}\\ \multicolumn{1}{c}{=\frac{{\partial }^{n}C(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_n\left(x_n\right);\theta )}{\partial F_1\left(x_1\right)\partial F_2\left(x_2\right)\cdots \partial F_n\left(x_n\right)}}\end{array}$$

(5)

The Copula function theory simplifies the construction of the joint distribution into two independent steps: the estimation of the marginal distribution and the selection of the Copula function. There are various kinds of Copula functions, such as Gaussian Copula function, t-Copula function, Gumbel Copula function, Clayton Copula function and Frank Copula function. Nelsen gives the function expression and Copula parameter range of these Copula functions. ²⁴ Because Gaussian Copula function only needs the marginal distribution and the correlation coefficient to determine the joint distribution uniquely, and can reflect the positive and negative correlation between variables (the range of correlation coefficient can reach (−1, 1)), ²¹ it is adopted to construct the prior joint distribution and the likelihood function in the Bayesian updating process of the aircraft lift device.

Reliability model of the high-lift device mechanism

Construct the wear model of the high-lift device mechanism

The high-lift device of civil airplanes is mainly made up of the pulley rack and sliding rail mechanism, the pinion and rack mechanism, the support arm and rocker arm mechanism, and so on. The pinion and rack mechanism is taken as the study object in this article, and its reliability is analyzed based on the wear mechanism. The sketch of the pinion and rack mechanism is shown in Figure 1.

Figure 1.

Pinion and rack mechanism.

Based on the Archard ²⁵ wear theory, the wear model of the pinion and rack mechanism is given by

$$W=K\frac{\mathrm{PS}\left(n_0\right)}{H}$$

(6)

where H is the material hardness, K is the wear coefficient, P is the normal load, $n_0$ is running times of the pinion and rack mechanism, and $S(n_0)$ is the total sliding distance, given by

$$S\left(n_0\right)=4n_0{L}_{0}S$$

(7)

in which the length of contact line $L_0$ is 19 mm, and S is the sliding distance of the pinion and rack mechanism running one time.

By looking up the metallic materials data handbook, K is $5\times {10}^{-9}$ , and H is 286 HB. P and S will be analyzed and calculated in next two sections.

Normal load P

Contact between the pinion and rack can be regarded as the contact of two gears whose axes are parallel from each other, but the diameters are different under a normal load P. The Hertz contact theory suggests that the normal load P is transferred through a point or line. Considering the deformation of the contact point after being loaded, the contact stress $\sigma$ can be calculated using the Hertz formula ²⁶

$$\sigma =\sqrt{\frac{P}{\pi }\frac{{\rho }_2-{\rho }_1}{{\rho }_1{\rho }_2(\frac{4G_1-E_1}{4{G_1}^2}+\frac{4G_2-E_2}{4{G_2}^2})L_0}}$$

(8)

where ${\rho }_1$ is the radius of contact point located in the pinion tooth profile; ${\rho }_2$ is the radius of contact point located in the rack tooth profile; $G_1$ and $G_2$ are the shear modulus of the pinion and rack, respectively; and $E_1$ and $E_2$ are the elastic modulus of the pinion and rack, respectively.

To simplify the calculation, assume that the main contact area of the pinion and rack locates in the pitch point in the meshing process of the pinion and rack. Therefore, the radii of pitch points of the pinion and rack can be acquired based on the involute equations and geometrical relationship of tooth profiles of the pinion and rack. The tooth profiles meshing schematic is shown in Figure 2. ²⁷

Figure 2.

Tooth profiles meshing.

The involute polar equations of the two contact tooth profiles are given by

$$r_i=r_{\mathrm{bi}}/\cos {\alpha }_0,i=1,2$$

(9)

where $r_1$ and $r_2$ are the radii of the pitch circle $O_1$ and $O_2$ , respectively; $r_{b1}$ and $r_{b2}$ are the radii of the base circle $O_1$ and $O_2$ , respectively; ${\alpha }_0$ is the pressure angle of the pitch circle $O_1$ and $O_2$ ; and ${\alpha }_0$ is 20°. From Figure 2, it is also seen that

$${\rho }_i=N_i{C}_0=r_i\sin {\alpha }_0,i=1,2$$

(10)

$r_1$ and $r_2$ are given by

$$r_i=\frac{1}{2}m{z}_i,i=1,2$$

(11)

in which m is the module, and m is 4.25; $z_1$ is the tooth number of the pinion, and $z_1$ is 13; $z_2$ is the tooth number of the rack, and $z_2$ is 245.

Substitute equations (10) and (11) into equation (8), we can obtain the normal load P based on the Hertz contact theory

$$P=\frac{1}{2}{\sigma }^2{L}_0\pi m\frac{z_1{z}_2}{z_2-z_1}(\frac{4G_1-E_1}{4{G_1}^2}+\frac{4G_2-E_2}{4{G_2}^2})\sin {\alpha }_0$$

(12)

Sliding distance S

The meshing process of the pinion and rack is the addendum of the pinion engaging-in the dedendum of the rack. Therefore, the sliding distance S can be regarded as the length of the tooth profile involute of the pinion approximately. Establish the involute equation of the pinion in the Cartesian coordinates as follows ²⁸

$$\{\begin{array}{c}x\left(\alpha \right)=r_{b1}\frac{\cos \left(\mathrm{inv}\left(\alpha \right)\right)}{\cos \left(\alpha \right)}\hfill \\ \multicolumn{1}{c}{y\left(\alpha \right)=r_{b1}\frac{\sin \left(\mathrm{inv}\left(\alpha \right)\right)}{\cos \left(\alpha \right)}}\end{array}$$

(13)

where $\alpha$ is the pressure angle, and $\mathrm{inv}(\alpha )$ is given by

$$\mathrm{inv}\left(\alpha \right)=\tan \left(\alpha \right)-\alpha$$

(14)

Substituting equation (14) into equation (13), we can obtain

$$\{\begin{array}{c}x\left(\alpha \right)=r_{b1}(\cos \left(\tan \left(\alpha \right)\right)+\sin \left(\tan \left(\alpha \right)\right)\tan \left(\alpha \right))\hfill \\ \multicolumn{1}{c}{y\left(\alpha \right)=r_{b1}(\sin \left(\tan \left(\alpha \right)\right)-\cos \left(\tan \left(\alpha \right)\right)\tan \left(\alpha \right))}\end{array}$$

(15)

By the integration of equation (15), we can obtain the sliding distance S (the length of involute)

$$S=\underset{0}{\overset{{\alpha }^{*}}{\int }}\sqrt{{\left(\frac{\mathrm{d}x\left(\alpha \right)}{\mathrm{d}\alpha }\right)}^2+{\left(\frac{\mathrm{d}y\left(\alpha \right)}{\mathrm{d}\alpha }\right)}^2}\mathrm{d}\alpha$$

(16)

where the pressure angle of the addendum circle of the pinion is given by

$${\alpha }^{*}=\arccos \left(\frac{r_{b1}}{r_{a1}}\right)=\arccos \left(\frac{z_1\cos {\alpha }_0}{z_1+2h_a^{*}}\right)$$

(17)

where addendum coefficient $h_a^{*}=1$ .

Substituting equation (17) into equation (16), we can see that the value of S is 9.5 mm. Then, substituting equation (7), equation (12), and the value of S into equation (6), we can obtain the wear model of the pinion and rack mechanism

$$W=a{n}_0{\sigma }^2(\frac{4G_1-E_1}{4G_1^2}+\frac{4G_2-E_2}{4G_1^2})$$

(18)

in which a is constant, and $a=1.3552\times {10}^{-6}$ .

Construct the wear reliability model of the high-lift device mechanism

In order to facilitate the calculation, the elastic modulus $E_2$ and shear modulus $G_2$ in equation (18) are assumed to be deterministic. We only consider the random uncertainty of the elastic modulus $E_1$ , the shear modulus $G_1$ , the contact stress $\sigma$ , and the permitted wear volume $W_0$ of the pinion and rack mechanism. Based on the stress-strength interference model, the pinion and rack mechanism fails when the actual wear volume $W(n_0,E_1,G_1,\sigma )$ exceeds the permitted wear volume $W_0$ . The wear reliability model of the pinion and rack mechanism is given by

$$g(n_0,E_1,G_1,\sigma ,W_0)=W_0-W(n_0,E_1,G_1,\sigma )$$

(19)

where $g(n_0,E_1,G_1,\sigma ,W_0)$ is the time-dependent performance function. Substituting equation (18) into equation (19), the wear reliability model is changed to

$$\begin{array}{c}\multicolumn{1}{c}{g(n_0,E_1,G_1,\sigma ,W_0)=W_0}\\ \multicolumn{1}{c}{-a{n}_0{\sigma }^2(\frac{4G_1-E_1}{4{G_1}^2}+\frac{4G_2-E_2}{4{G_2}^2})}\end{array}$$

(20)

Actually, the elastic modulus $E_1$ and shear modulus $G_1$ are not independent from each other. The function relation between them is given by

$$E_1=2(1+{\nu }_1)G_1$$

(21)

in which ${\nu }_1$ is Poisson’s ratio.

The correlation between $E_1$ and $G_1$ can be described by the Pearson correlation coefficient. It is necessary to study the impact of the correlation on the reliability analysis results. Otherwise, neglecting the correlation of random variables may cause inaccurate results.

Reliability analysis of the high-lift device mechanism

Update cases

Case 1: update both E₁ and G₁ separately

In this case, both $E_1$ and $G_1$ are updated, but they are considered as independent from each other. In case 2, both $E_1$ and $G_1$ are updated simultaneously, and the correlation between them is considered. By comparing the two cases, we can see the impact of the correlation on updating and reliability analysis results.

Assume that both of the prior probability distributions of $E_1$ and $G_1$ follow the uniform distribution. They are selected as respectively: $E_1~\mathrm{U}\left(180000,200000\right)$ and $G_1~\mathrm{U}(73000,89000)$ . Based on the distribution information, the mean value, standard deviation, and coefficient of variation of $E_1$ are, respectively, 200,000 MPa, 11,547 MPa, and 0.058; the mean value, standard deviation, and coefficient of variation of $G_1$ are, respectively, 81,000 MPa, 4619 MPa, and 0.057. The prior PDFs of $E_1$ and $G_1$ are given by, respectively

$$f_1\left(E_1\right)=\{\begin{array}{cc}\frac{1}{40000},\hfill & 180,000\le E_1\le 220,000\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{otherwise}\hfill \end{array}$$

(22)

$$f_2\left(G_1\right)=\{\begin{array}{cc}\frac{1}{16,000},\hfill & 73,000\le G_1\le 89,000\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{otherwise}\hfill \end{array}$$

(23)

Besides the prior PDF, we also need the likelihood function to calculate the posterior PDF. The experimental data are obtained from the research project with the First Aircraft Institute: Reliability Analysis of High-lift Device Mechanism. In this research project, the measurement experiment about the elastic modulus, Poisson’s ratio, and shear elasticity of materials used for making the high-lift device mechanism is carried out.

Based on the above experimental data, the mean value ${\mu }_{L1}$ and standard deviation ${\sigma }_{L1}$ of elasticity modulus are, respectively, 205,676 and 3742 MPa; the mean value ${\mu }_{L2}$ and standard deviation ${\sigma }_{L2}$ of shear modulus are, respectively, 80,944 and 1449 MPa. Also, we can infer from the experimental data that the likelihood functions of $E_1$ and $G_1$ are both approximately follow the normal distribution, given by

$$L_1\left(E_1\right)=\frac{1}{2\pi {\sigma }_{L1}}\exp (-\frac{{(E_1-{\mu }_{L1})}^2}{2{\sigma }_{L1}^2})$$

(24)

$$L_2\left(G_1\right)=\frac{1}{2\pi {\sigma }_{L2}}\exp (-\frac{{(G_1-{\mu }_{L2})}^2}{2{\sigma }_{L2}^2})$$

(25)

According to equations (1) and (2), the posterior PDFs of $E_1$ and $G_1$ can be obtained. Then, the MCMC method is used to generate the posterior samples of $E_1$ and $G_1$ . Figures 3 and 4 show the sampling results of $E_1$ and $G_1$ .

Figure 3.

Trace of iteration in case 1: (a) trace of iteration of E₁ and (b) trace of iteration of G₁.

Figure 4.

Histogram and fitted PDFs of the generated posterior samples in case 1: (a) PDF of E₁ and (b) PDF of G₁.

Figure 3(a) represents the traces of 10,000 samples of $E_1$ . The initial samples are given up, since they do not converge. From Figure 3(a), we can see that the value in which the samples are given up is about 500 in this case. The mean value of the posterior distribution of $E_1$ converges to the neighborhood of 205,000 MPa with the increase in the iterations. Figure 3(b) represents traces of 10,000 samples of $G_1$ , and the value where the samples are given up is about 1000. The mean value of the posterior distribution of $G_1$ is about 80,800 MPa as shown in Figure 3(b).

Figure 4 shows the estimated PDFs of $E_1$ and $G_1$ with the samples. We can see that though the prior distributions of $E_1$ and $G_1$ both follow the uniform distribution, their posterior distributions both follow the normal distribution approximately.

Based on the above posterior samples, the mean value, standard deviation, and coefficient of variation of $E_1$ are, respectively, 205,589 MPa, 4898 MPa and 0.0238; the mean value, standard deviation, and coefficient of variation of $G_1$ are, respectively, 80,988 MPa, 2493 MPa and 0.0308 in this case.

Comparing the posterior with the prior, it is found that:

1. The distributions of $E_1$ and $G_1$ both change from the assumed uniform distribution to the normal distribution, as shown in Figure 4;

2. The mean values of $E_1$ increase, but the mean values of $G_1$ decrease slightly;

3. The standard deviations and coefficients of variation of $E_1$ and $G_1$ decrease obviously.

We can see easily from the above point 3 that the integration of experimental data effectively reduces the epistemic uncertainty of $E_1$ and $G_1$ . Because our knowledge about $E_1$ and $G_1$ is improved after the Bayesian updating is finished. The differences between the prior and posterior probability distributions of $E_1$ and $G_1$ are significant.

Case 2: update E₁ and G₁ simultaneously, considering the correlation

Different from case 1, $E_1$ and $G_1$ in this case are updated simultaneously and the correlation between them is taken into account during the updating process:

• 1. Construct the prior joint probability distribution using the Gaussian Copula function.

Based on the uniform distribution assumption in case 1, the prior marginal CDFs of $E_1$ and $G_1$ are given by, respectively

$$F_1\left(E_1\right)=\{\begin{array}{cc}0\hfill & E_1<180,000\hfill \\ \multicolumn{1}{c}{\frac{E_1-180,000}{40,000}}& 180,000\le E_1<220,000\hfill \\ \multicolumn{1}{c}{1}& E_1\ge 220,000\hfill \end{array}$$

(26)

$$F_2\left(G_1\right)=\{\begin{array}{cc}0\hfill & G_1<73,000\hfill \\ \multicolumn{1}{c}{\frac{G_1-73,000}{16,000}}& 73,000\le G_1<89,000\hfill \\ \multicolumn{1}{c}{1}& G_1\ge 89,000\hfill \end{array}$$

(27)

The Gaussian Copula density function of $E_1$ and $G_1$ is given as ²⁴

$$\begin{array}{c}\multicolumn{1}{c}{c(F_1\left(E_1\right),F_2\left(G_1\right);\theta )=\frac{1}{\sqrt{1-{\theta }^2}}}\\ \multicolumn{1}{c}{\exp (-\frac{y_1^2{\theta }^2-2\theta y_1{y}_2+y_2^2{\theta }^2}{2(1-{\theta }^2)})}\end{array}$$

(28)

in which $y_1={\mathrm{\Phi }}^{-1}(F_1(E_1))$ and $y_2={\mathrm{\Phi }}^{-1}(F_2(G_1))$ are both the standard normal random variables.

To acquire the prior joint PDF of $E_1$ and $G_1$ , we also should know the Copula parameter $\theta$ . Assume that the symbol ${\rho }_P$ represents the Pearson correlation coefficient of $E_1$ and $G_1$ in the prior probability distribution. The function relationship between the Pearson correlation coefficient ${\rho }_P$ and the Copula parameter $\theta$ is given by ¹⁵

$$\begin{array}{l}{\rho }_P=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}\left(\frac{E_1-{\mu }_1}{{\sigma }_1}\right)\left(\frac{G_1-{\mu }_2}{{\sigma }_2}\right)\frac{f_1\left(E_1\right)f_2\left(G_1\right)}{\sqrt{1-{\theta }^2}}\\ \exp \left(-\frac{y_1^2{\theta }^2-2\theta y_1{y}_2+y_2^2{\theta }^2}{2\left(1-{\theta }^2\right)}\right)\mathrm{d}{E}_1\mathrm{d}{G}_1\end{array}$$

(29)

in which ${\mu }_1$ and ${\mu }_2$ represent the mean value of $E_1$ and $G_1$ , respectively; ${\sigma }_1$ and ${\sigma }_2$ represent the standard deviation of $E_1$ and $G_1$ , respectively; $f_1(E_1)$ and $f_2(G_1)$ represent the prior marginal PDF of $E_1$ and $G_1$ , respectively. A MATLAB program is written to calculate $\theta$ based on equation (29) in this article. Assuming that the Pearson correlation coefficient ${\rho }_P$ between $E_1$ and $G_1$ in the prior joint PDF is 0.5, we can obtain the value of Copula parameter $\theta$ by the MATLAB program, and the value of $\theta$ is 0.52. Then, substituting equations (22), (23), (26)–(28), and the value of $\theta$ into equation (4), the prior joint PDF of $E_1$ and $G_1$ can be determined.

• 2. Construct the likelihood function using the Gaussian Copula function.

We also use the Gaussian Copula density function to construct the likelihood function of $E_1$ and $G_1$ . Based on the experimental data, the Pearson correlation coefficient ${\rho }_{\mathrm{PL}}$ between the elasticity modulus and shear modulus is 0.94. Therefore, the value of the Copula parameter ${\theta }_L$ is 0.96 by the MATLAB program. The likelihood function of $E_1$ and $G_1$ considering their correlation is given as

$$L(E_1,G_1)=c_L(F_{L1}\left(E_1\right),F_{L2}\left(G_1\right);{\theta }_L)L_1\left(E_1\right)L_2\left(G_1\right)$$

(30)

in which $F_{L1}(E_1)$ , $F_{L2}(G_1)$ , and $c_L(·)$ are given by

$$\begin{array}{c}\multicolumn{1}{c}{F_{L1}\left(E_1\right)=\underset{-\infty }{\overset{E_1}{\int }}{L}_1\left(E_1\right)\mathrm{d}{E}_1}\\ \multicolumn{1}{c}{=\underset{-\infty }{\overset{E_1}{\int }}\frac{1}{2\pi {\sigma }_{L1}}\exp (-\frac{{(E_1-{\mu }_{L1})}^2}{2{\sigma }_{L1}^2})\mathrm{d}{E}_1}\end{array}$$

(31)

$$\begin{array}{c}\multicolumn{1}{c}{F_{L2}\left(G_1\right)=\underset{-\infty }{\overset{G_1}{\int }}{L}_2\left(G_1\right)\mathrm{d}{G}_1}\\ \multicolumn{1}{c}{=\underset{-\infty }{\overset{G_1}{\int }}\frac{1}{2\pi {\sigma }_{L2}}\exp (-\frac{{(G_1-{\mu }_{L2})}^2}{2{\sigma }_{L2}^2})\mathrm{d}{G}_1}\end{array}$$

(32)

$$\begin{array}{c}\multicolumn{1}{c}{c_L(F_{L1}\left(E_1\right),F_{L2}\left(G_1\right);{\theta }_L)=\frac{1}{\sqrt{1-{\theta }_L^2}}}\\ \multicolumn{1}{c}{\exp (-\frac{y_{L1}^2{\theta }_L^2-2{\theta }_L{y}_{L1}{y}_{L2}+y_{L2}^2{\theta }_L^2}{2(1-{\theta }_L^2)})}\end{array}$$

(33)

where $y_{L1}={\mathrm{\Phi }}^{-1}(F_{L1}(E_1))$ and $y_{L2}={\mathrm{\Phi }}^{-1}(F_{L2}(G_1))$ are both the standard normal random variables.

Then, substituting equations (31)–(33) and the value of ${\theta }_L$ into equation (30), the likelihood function of $E_1$ and $G_1$ considering correlation can be determined.

• 3. Update the prior joint PDF.

With the prior joint PDF and the likelihood function, the posterior samples of $E_1$ and $G_1$ can be obtained based on the MCMC method, as plotted in Figures 5 and 6.

Figure 5.

Trace of iteration in case 2: (a) trace of iteration of E₁ and (b) trace of iteration of G₁.

Figure 6.

Histogram and fitted PDFs of the generated posterior samples in case 2: (a) PDF of E₁ and (b) PDF of G₁.

Figure 5(a) shows the traces of 10,000 samples of $E_1$ . From Figure 5(a), we can see that the value in which the samples are given up is about 1000 in this case. The mean value of the posterior distribution of $E_1$ converges to the neighborhood of 205,000 MPa with the increase in the iterations, which is similar to the observation in case 1. Figure 5(b) shows the traces of 10,000 samples of $G_1$ , and the value where the samples are given up is also about 1000. The mean value of the posterior distribution of $G_1$ is about 80,900 MPa.

Figure 6 shows the estimated posterior marginal PDFs of $E_1$ and $G_1$ . Different from case 1, the posterior marginal probability distributions both follow the lognormal distribution.

In this case, according to the posterior samples, the mean value, standard deviation, and coefficient of variation of $E_1$ in the posterior probability distribution are, respectively, 205,665 MPa, 3814 MPa and 0.0185; the mean value, standard deviation, and coefficient of variation of $G_1$ in the posterior probability distribution are, respectively, 80,956 MPa, 1497 MPa and 0.0185; and the Pearson correlation coefficient between $E_1$ and $G_1$ is 0.85. By comparison, it is found that:

1. The posterior marginal probability distributions of $E_1$ and $G_1$ both change to the lognormal distribution; the observations are different from those in the prior probability distributions (uniform distribution) and those in case 1 (normal distribution);

2. Similar to the observations in case 1, the mean value of $E_1$ increases, and the standard deviation and coefficient of variation of $E_1$ decrease obviously in this case compared to the prior probability distribution of $E_1$ ; the reason is (1) Bayesian updating makes the mean value and standard deviation of $E_1$ come close to those associated with the experimental data; (2) the mean value of the experimental data is larger than the mean value of the prior probability distribution of $E_1$ ; however, the standard deviation of the experimental data is smaller than that of the prior probability distribution of $E_1$ ;

3. Similar to the observations in case 1, the mean values of $G_1$ decrease slightly, and the standard deviation and coefficient of variation of $G_1$ decrease obviously in this case compared to the prior probability distribution of $G_1$ ; the reason for this is (a) Bayesian updating makes the mean value and standard deviation of $G_1$ come close to those associated with the experimental data; (b) the mean value and standard deviation of the experimental data are both smaller than those of the prior probability distribution of $G_1$ ;

4. Compared with case 1, the differences in the standard deviations and coefficients of variation between the prior and posterior in this case are more significant. The reason is that the two random variables $E_1$ and $G_1$ are actually correlated. Therefore, the independence assumption in case 1 is not reasonable, and the obtained posterior marginal PDFs in this case are closer to the true results than those in case 1;

5. The correlation between $E_1$ and $G_1$ increases after incorporation of the experimental data;

6. The posterior joint PDF is also different from the prior joint PDF, as shown in Figure 7.

Figure 7.

Prior and posterior joint PDFs of E₁ and G₁: (a) prior joint PDF of E₁ and G₁ and (b) posterior joint PDF of E₁ and G₁.

Reliability analysis

According to equation (20), the performance function $g(n_0,E_1,G_1,\sigma ,W_0)$ changes when the running times of the high-lift device change. This is essentially a time-dependent reliability problem. If $N_0$ is the lifetime of interest, then the reliability within the lifetime $[0,N_0]$ can be described as

$$R(0,N_0)=\mathrm{Prob}(g(n_0,E_1,G_1,\sigma ,W_0)>0,\exists n_0\in [0,N_0])$$

(34)

Before calculating $R(0,N_0)$ , let us first define the so-called instantaneous reliability $R(n_0)$ at $n_0$ by

$$R\left(n_0\right)=\mathrm{Prob}(g(n_0,E_1,G_1,\sigma ,W_0)>0)$$

(35)

The instantaneous reliability $R(n_0)$ is calculated by considering running times as a fixed parameter in equation (35). Note that this probability does not take into account what happened before $n_0$ . It is independent of the previous states, that is, it is not a conditional probability. It will be used in the following.

It is extremely difficult to analyze $R(0,N_0)$ directly because it requires numerical analysis of multidimensional integration. Based on the equivalent extreme value approach, ²⁹ the pinion and rack mechanism is reliable as long as $g_{min}>0$ , where $g_{min}$ is the global minimum value of the performance function $g(n_0,E_1,G_1,\sigma ,W_0)$ on $[0,N_0]$ is given by

$$g_{min}=\underset{n_0\in [0,N_0]}{\arg min}\left(g(n_0,E_1,G_1,\sigma ,W_0)\right)$$

(36)

Then, $R(0,N_0)$ can be rewritten as

$$R(0,N_0)=\mathrm{Prob}(g_{min}>0)$$

(37)

According to equations (36) and (37), the above time-dependent reliability problem can be converted to a time-independent one. The key to realize this conversion is to calculate the global minimum value $g_{min}$ . Note that the actual wear volume $W(n_0,E_1,G_1,\sigma )$ increases monotonically over running times $n_0$ , so the global minimum value of the performance function $g(n_0,E_1,G_1,\sigma ,W_0)$ will generally occur at the running times–interval boundary. Therefore, the reliability calculation needs to be carried out only at the running times–interval boundary, which can guarantee the reliability being satisfied over $[0,N_0]$ .

In summary, the calculation of $R(0,N_0)$ as defined in equation (34) reduces to that of $R(N_0)$ as defined in equation (35). Equation (34) is changed to

$$\begin{array}{cc}\multicolumn{1}{c}{R(0,N_0)}& =\mathrm{Prob}(g(n_0,E_1,G_1,\sigma ,W_0)>0,\exists n_0\in [0,N_0])\hfill \\ \multicolumn{1}{c}{}& =\mathrm{Prob}(g(N_0,E_1,G_1,\sigma ,W_0)>0)\hfill \\ \multicolumn{1}{c}{}& =R\left(N_0\right)\hfill \end{array}$$

(38)

Thus, the existing reliability analysis approach, such as the first-order second moment (FOSM) method, the second-order second moment (SOSM) method, and Monte Carlo (MC) method, can be conveniently used.

Assume that the contact stress $\sigma$ and the permitted wear volume $W_0$ both follow the normal distribution, and the probability distribution parameters are given in Table 1. Then, we can obtain the reliability analysis results of the pinion and rack mechanism based on equations (35) and (38). Finally, we can obtain the generalized reliability indices. The generalized reliability index is defined by formula ³⁰

$$G\beta =-{\mathrm{\Phi }}^{-1}(1-R(0,N_0))$$

(39)

Table 1.

Probability distributions of W₀ and σ.

Variablename	Mean value	Standarddeviation	Distributiontype
σ	480 MPa	48 MPa	Normal
W 0	0.03	0.003	Normal

The generalized reliability indices of the pinion and rack mechanism associated with cases 1 and 2 are plotted in Figure 8. It is seen that:

1. The prior generalized reliability index (Gbeta01 curve, before updating and without considering correlation) is the lowest;

2. The generalized reliability index associated with updating the $E_1$ and $G_1$ simultaneously (Gbeta03 curve, considering correlation) is the highest;

3. Because of the integration of the partial data (containing no correlation information), the generalized reliability index associated with updating both $E_1$ and $G_1$ separately (Gbeta02 curve, without considering correlation) is lower than the that associated with updating the $E_1$ and $G_1$ simultaneously; the gap between the two generalized reliability indices first widens as the $N_0$ increases. When $N_0=5000$ , the gap reaches its maximum. The maximum is about 10.4%. Then, the gap gradually narrows; when $N_0=8000$ , it changes to about 7.7%. The gap is summarized in Table 2.

Figure 8.

The generalized reliability indices before and after updating.

Table 2.

The gap of generalized reliability indices associated with case 1 and case 2.

N 0	1000	2000	3000	4000	5000	6000	7000	8000
Gap of generalized reliability indices (%)	3.4	5.4	8.1	9.8	10.4	10.1	9.1	7.7

According to the above analysis, we can see that the differences between the posterior probability distribution considering correlation and that without considering correlation is obvious. Because $E_1$ and $G_1$ are actually correlated, the posterior probability distribution considering correlation is closer to the true result. Only when we get accurate probability distributions, then the credible reliability analysis results can be obtained. Otherwise, we may obtain inaccurate analysis results as shown in Table 2.

Conclusion

1. Under the condition of the incomplete probability information, this article presents an approach for improving the reliability analysis accuracy of the aircraft high-lift by introducing the Copula function theory into Bayesian updating. The prior joint PDF and likelihood function of correlated random variables are constructed by the Copula function. The posterior joint PDF is obtained by Bayesian updating.

2. Taking the pinion and rack mechanism of the aircraft high-lift device for an example, two cases associated with updating (a) two random variables separately and (b) two correlated random variables simultaneously are investigated. The analysis of two cases shows that the reliability results based on the proposed approach are more accurate and more coincident with the factual situation. This approach can also be applied to other products’ reliability analysis and extend the application scope of Bayesian updating effectively.