Time-dependent system reliability is computed as the probability that the responses of a system do not exceed prescribed failure thresholds over a time duration of interest. In this work, an efficient time-dependent reliability analysis method is proposed for systems with bivariate responses which are general functions of random variables and stochastic processes. Analytical expressions are derived first for the single and joint upcrossing rates based on the first-order reliability method (FORM). Time-dependent system failure probability is then estimated with the computed single and joint upcrossing rates. The method can efficiently and accurately estimate different types of upcrossing rates for the systems with bivariate responses when FORM is applicable. In addition, the developed method is applicable to general problems with random variables, stationary, and nonstationary stochastic processes. As the general system reliability can be approximated with the results from reliability analyses for individual responses and bivariate responses, the proposed method can be extended to reliability analysis of general systems with more than two responses. Three examples, including a parallel system, a series system, and a hydrokinetic turbine blade application, are used to demonstrate the effectiveness of the proposed method.

## Introduction

Reliability is the ability that a component or system performs its intended function in routine circumstances for a given period of time [1]. For a system, we are interested in not only the reliability of its components, but also the system reliability. For many applications, the system reliability is more critical, because it can be used to estimate the lifecycle cost, determine maintenance policies, and perform optimization [2–4]. The system reliability is determined by the reliability of its components and also the relationship between the system and its components. By component, in this work, we mean a failure mode. A physical component itself can also be regarded as a system, because it may have multiple failure modes.

System reliability analysis is much more difficult than component reliability analysis. The latter needs to compute the probability of a single event, while the former needs to compute the joint probability of multiple and possibly dependent events. Many progresses have been made in system reliability analysis. For example, Ditlevsen [5] approximated the system reliability using a bounding formulas. Song and Kang [6] developed a matrix-based system reliability method, which can calculate the system reliability and system parameter sensitivities by a convenient matrix-based framework. Nguyen et al. [7] later developed a reliability-based system design optimization method by using the matrix-based system reliability method. Dey and Mahadevan [8] and Ambartzumian et al. [9] proposed a system reliability method using a standard normal multivariate cumulative distribution function (CDF); by employing the Morgan's laws [10], the method expresses the system probability of failure as the intersection of a set of unions of subsystems. Youn and Wang [11] suggested a complementary intersection method (CIM) for system reliability analysis by expressing the system failure event as complementary intersection events. Based on the CIM method, a generalized CIM method is developed later in Ref. [12]. More system reliability analysis methods have also been reported in Refs. [1] and [13].

To predict the system reliability, we need at first to calculate component reliability. This can be achieved by using the common reliability methods, such as the first-order reliability method (FORM), second-order reliability method, and importance sampling [1,14–19]. The next step is to calculate joint probabilities associated with the failure events of all the components. If there are only two failure modes, the system reliability can be computed using a bivariate distribution. One example was reported in Ref. [9]. When there are many components involved, it is difficult and computationally expensive to obtain the joint probabilities. For this case, the system reliability can be approximated within its lower and upper bounds by just using the component reliabilities and joint probabilities for up to two failure modes (bivariate responses) [5]; the reliability bound is normally narrow for many engineering problems.

Although reliability is defined for a period of time and is also a function of time, most of the aforementioned reliability methods are for time-invariant reliability that does not change over time. For time-invariant reliability, the associated function of a response variable, called a limit-state function, is not a function of time. The reliability is therefore time independent. In many engineering applications, however, the limit-state function changes over time, either because time appears explicitly in the function or because stochastic processes are part of the input variables, or both. Examples include function generator mechanisms [20,21], bridges under stochastic loading [22,23], hydrokinetic turbine system subjected to wave or river flow loading [24,25], and vehicles running on stochastic road surfaces [4].

Time-dependent reliability analysis is much more challenging than its time-independent counterpart. The most common time-dependent reliability method is the Rice formula [26,27] developed in 1944 and is still widely used nowadays. There are still many developments in time-dependent reliability in recent years. For instance, for component reliability problems, Yang and Shlnozuka [28,29] developed a time-dependent reliability analysis method using the Markov stochastic process. Mourelatos and coworkers [30] employed the time-series modeling and importance sampling method to approximate the time-dependent reliability. Andrieu-Renaud et al. [31] proposed a PHI2 method for the time-dependent component reliability analysis of general problems with both random variables and stochastic processes. By using the Rice's formula [26,27], Du and Hu [24] developed a time-dependent reliability model for hydrokinetic turbine blades. Singh et al. [4] proposed the concept of composite limit-state function for time-dependent reliability analysis for a special group of problems. Du and Hu [32] presented a joint upcrossing method based on the work of Madsen and Krenk [33]. Wang and Wang [34,35] developed a nested extreme value response surface method by building a surrogate model for the extreme value response. They also proposed a double-loop surrogate modeling approach for time-dependent component reliability analysis in Ref. [36]. Jiang et al. [37,38] studied the time-dependent reliability analysis of general engineering systems based on stochastic process discretization.

Studies on time-dependent system reliability have also been reported. For example, a method was developed to estimate the joint first-passage probability of failure for systems under stochastic excitation [39]. An analytical approach was proposed for linear dynamical systems in higher dimensions by calculating the conditional upcrossing rate on the surfaces of failure boundaries [40]. An approximation method was reported to approximate the conditional first passage probability of systems under modulated white noise excitation [41]. By combining Monte Carlo simulations (MCS) with the asymptotic extreme value theory, Radhika et al. [42] proposed a reliability analysis method for nonlinear vibrating systems. Gupta and Manohar [43] proposed a multivariate extreme value distribution approach to approximate the extreme value distributions of a vector of stationary Gaussian random processes. The developed multivariate extreme value distribution method was then applied to reliability analysis of randomly vibrating structures subjected to jointly stationary Gaussian excitations [44,45]. Some of the above methods have been verified to have good accuracy for systems subjected to multiple Gaussian stationary stochastic processes. They, however, cannot be directly applied to general problems where the input variables of a limit-state function contain time, random variables, and nonstationary stochastic processes. Hagen and Tvedt [46] have investigated the reliability analysis of parallel system under vector of stochastic processes. Their method, however, needs to solve multivariate normal distribution, and new random variables are introduced into the model to compute the outcrossing rate. Although MCS is capable of handling the general problems, the required computational cost is prohibitive for engineering applications.

In this work, we propose a new time-dependent system reliability analysis method for bivariate responses that are nonlinear functions of time, random variables, and stochastic processes. The stochastic process can be stationary or nonstationary. The new method is an extension of the work in Ref. [39] and is based on the FORM and the Rice's formula [26,27]. The major development is the derivations of bivariate joint upcrossing rates, which can be used for estimating time-dependent system reliability for bivariate responses. Since the bivariate joint probabilities are the basis for general system reliability analysis when the reliability bound method is used, the proposed method can also be applied to general time-dependent system reliability analysis for general systems with more than two components.

In Sec. 2, we give the background of time-dependent reliability. In Sec. 3, we first introduce the upcrossing rate method; we then derive necessary equations for the bivariate joint upcrossing rates. The numerical procedure is summarized in Sec. 4, followed by three examples in Sec. 5. Conclusions are presented in Sec. 6.

## Time-Dependent System Reliability

Recall that in this work, a component corresponds to a failure mode. Suppose there are *r* failure modes or *r* components, for component *i,* where $i=1,2,\u2026,r$, let its limit-state function be $Gi(t)=gi(X,\u2009Y(t),\u2009t)$, where $X=[X1,\u2009X2,\u2009\cdots ,\u2009Xn]$ is a vector of random variables, $Y(t)=[Y1(t),\u2009Y2(t),\u2009\cdots ,\u2009Ym(t)]$ is a vector of stochastic processes, $Gi(t)$ is the response variable, and *t* stands for time. Functions, $Gi(t)=gi(X,\u2009Y(t),\u2009t)$, $i=1,\u20092,\u2009\cdots ,\u2009r$, are usually computer simulation models, such as those of finite element analysis (FEA) or computational fluid dynamics. They are nonlinear functions of random variables $X$ and stochastic processes $Y(t)$. As a result, the response variables $Gi(t)$, $i=1,\u20092,\u2009\cdots ,\u2009r$, are stochastic processes with variation and correlations over time. For a special case where $Gi(t)=gi(X,\u2009t)$, although no stochastic processes $Y(t)$ appear in the input, the response variable $Gi(t)$ is still a stochastic process due to the time-varying statistical properties and correlations over time. In addition, even if $X$ and $Y(t)$ are mutually independent from each other, the stochastic responses $Gi(t)$, $i=1,\u20092,\u2009\cdots ,\u2009r$, are still correlated multivariate stochastic processes due to the sharing random variables or stochastic processes in the limit-state functions.

*i*over the time interval $[t0,\u2009ts]$ is defined by

*r*. Evaluating a joint probability with a high-order is extremely difficult. To make the system reliability easier, Ditlevsen [5] proposed abound formulas for a series system as follows:

In the above equations, $pf,i(t0,\u2009ts)$ is the component probability of failure. Components probabilities of failure are sorted in a decreasing order, and $pf,1(t0,\u2009ts)$ is, therefore, the maximum component probability of failure [1]. $pf,ij(t0,\u2009ts)$ is the joint probability of failure of components *i* and *j*. It is the probability that both components fail over $[t0,\u2009ts]$. The above reliability bounds (i.e., Eqs. (7)–(9)) are for series systems. For parallel systems, the failure probability can be expressed as a function of series system reliability using De Morgan's law [10]. For mixed system with combined series and parallel structures, the generalized CIM method presented in Ref. [12] can be employed to decompose the system failure event into multiple mutually exclusive failure paths. Based on the decomposition, the reliability bounds for series systems can be used to provide bounds for the system reliability.

Equation (9) indicates that the component probability of failure $pf,i(t0,\u2009ts)$ and bivariate probability of failure $pf,ij(t0,\u2009ts)$ are the bases for the system probability of failure. As reviewed previously, many time-dependent reliability methods are available for $pf,i(t0,\u2009ts)$. In this work, we focus on developing a new method for the bivariate probability of failure $pf,ij(t0,\u2009ts)$. A straightforward way of evaluating $pf,ij(t0,\u2009ts)$ is to directly generate samples of $X$ and $Y(t)$ using Monte Carlo simulation method. After that, the samples of $X$ and $Y(t)$ are plugged into the limit-state functions to get samples of responses and then obtain the system failure probability. Estimating the failure probability in this way is straightforward and accurate. However, as discussed previously, the limit-state functions may be expensive simulation models, such as FEA and computational fluid dynamics models. The required computational effort for this kind of straightforward simulation method is prohibitive. We discuss how to overcome this challenge in Sec. 3. The idea is to reduce the number of times of calling the limit-state function while maintaining the accuracy of reliability analysis.

## Time-Dependent System Reliability for Bivariate Responses

*i*and

*j*. In Secs. 3.1–3.3, we first discuss the method for the time-dependent component reliability. We then derive equations for time-dependent joint probability, $pf,\u2009i\u2009\u222a\u2009j(t0,\u2009ts)$.

### Time-Dependent Component Reliability Analysis.

In this work, we employ the upcrossing rate method [31] to evaluate the time-dependent component probability of failure. It is the most commonly used method for time-dependent reliability analysis.

#### Upcrossing Rate Method for Time-Dependent Component Reliability Analysis.

*k*at time instant

*t*. Figure 1 shows upcrossing events of component

*k*. An upcrossing event happens when the response variable

*G*passes the threshold

_{k}*e*at time instant

_{k}*t*from the safe region $Gk(t)<ek$ to the failure region $Gk(t+\Delta t)>ek$, where $\Delta t$ is an infinitesimally small time interval. The upcrossing rate $vk+(t)$ is defined by the following limit:

Equation (14) is derived based on the assumption that all the upcrossings over $[t0,\u2009ts]$ are independent. Knowing $vk+(t)$, we can easily obtain the component probability of failure $pf,k(t0,\u2009ts)$ using Eq. (14).

In addition to using the upcrossing rate method, we also use FORM, which linearizes the limit-state function of components *i* and *j* at the so-called most probable point (MPP). Next, we first discuss the linearization and then discuss the estimation of upcrossing rate $vk+(t)$. The linearization later is also used in the derivation of $pf,\u2009i\u2009\u222a\u2009j(t0,\u2009ts)$.

#### Transformation of Limit-State Functions.

*k,*$Gk(t)=gk(X,\u2009Y(t),\u2009t),\u2009\u2009k=i\u2002or\u2002j$, we first transform random variables $X$ and stochastic processes $Y(t)$ into standard normal random variables $U(t)=(UX,\u2009UY(t))$. Then, the limit-state function becomes

where $T(\u22c5)$ stands for the transforming operator. The transformation can be found in Ref. [47].

$\beta k(t),\u2009k=i\u2009or\u2009j$ is called the Hasofer-Lind reliability index.

The equations for the time-dependent component probabilities of failure, $pf,\u2009k(t0,\u2009ts),\u2009k=i\u2009or\u2009j$, are already available. We reviewed these equations in Appendix A. In Sec. 3.2, we investigate the method for the approximation of the bivariate probability $pf,\u2009i\u2009\u222a\u2009j(t0,\u2009ts)$. Note that since the FORM method is employed to linearize the limit-state function at the MPP, the accuracy of the failure probability estimate may be affected by the accuracy of FORM. The developed method in Secs. 3.2 and 3.3 is therefore applicable to problems where FORM is accurate.

### Time-Dependent Joint Probability $pf,\u2009i\u2009\u222a\u2009j(t0,\u2009ts)$

#### Outcrossing Rate Method for Time-Dependent Joint Probability Analysis.

*i*and

*j*at time instant

*t*. An outcrossing event occurs when the system outcrosses its bounds at time instant

*t*from the safe region to the failure region. Figure 2 shows three representative outcrossing events of the series system. For the outcrossing events, both components

*i*and

*j*are in the safe region at time instants

*t*= 1, 2, and 3. The system then outcrosses into the failure region as a result of the upcrossing of

_{m}, m*G*, or upcorring of

_{i}*G*, or both the upcrossings of

_{j}*G*and

_{i}*G*at the following time instants,

_{j}*t*△

_{m}+*t, m*= 1,3, and 2. Given in mathematical form, the outcrossing rate $vi\u2009\u222a\u2009j+(t)$ is given by the following limit:

where $\Delta t$ is an infinitesimally small time interval.

$pij+\u2212(t)$ is the probability that $Gi(t)$ upcrosses its barrier $ei$, while $Gj(t)$ remains below its barrier $ej$ at $t$, $pij\u2212+(t)$ is the probability that $Gj(t)$ upcrosses its barrier $ej$, while $Gi(t)$ remains below its barrier $ei$ at $t$, and $pij++(t)$ is the probability that both $Gi(t)$ and $Gj(t)$ upcross their barriers at $t$.

Equations for $vij+\u2212(t)$ are available for special limit-state functions with two stationary Gaussian vector processes [39]. In Secs. 3.2.2–3.2.4, we derive equations for $vij+\u2212(t)$ and other two joint upcrossing rates for general limit-state functions. The derivations are based on the approximation discussed in Sec. 3.1.2.

#### $vij+\u2212(t)$.

where $fLiLjL\u02d9i(\u22c5,\u2009\u22c5,\u2009\u22c5)$ is the joint probability density function (PDF) of $Li(t)$, $Lj(t)$, and $L\u02d9i(t)$.

As there is no closed-form available for $vij+\u2212(t)$, we perform some transformations for Eq. (37) before we derive necessary equations for it. The transformation is given in Appendix B.

The above equations indicate that $\mu Lj|Li=\beta i$, $\sigma Lj|Li=\beta i$, $\mu L\u02d9i|Li=\beta i(t),\u2009Lj=lj$, and $\sigma L\u02d9i|Li=\beta i(t),\u2009Lj=lj$ are required to solve for $vij+\u2212(t)$, for which we must to obtain the mean and covariance of $L\u02d9i(t)$ and $L=[Li(t),\u2009Lj(t)]$.

#### $vij\u2212+(t)$.

#### $vij++(t)$.

*t*is omitted in Eq. (55) for brevity; for example, $zi$ stands for $zi(t)$ now. In Appendix C, we have demonstrated that $vij++(t)$ tends to be zero when $\Delta t$ becomes infinitely small. Based on the demonstration, we, therefore, conclude that $vij++(t)=0$.

#### $Rij(t0)$.

*i*and

*j*are found at $t0$, $Rij(t0)$ is calculated by [8]

### System Reliability Analysis.

With the availability of the outcrossing rate $vi\u2009\u222a\u2009j+(t)$, we now summarize the system reliability analysis method for bivariate responses.

In above equations, $Rij(t0)$ is computed using Eqs. (57) and (58), $vi\u2009\u222a\u2009j+(t)$ is estimated using Eq. (56), and $vi+(t)$ is computed using Eq. (16). In addition, $vij+\u2212(t)$ and $vij\u2212+(t)$ are required to compute $vi\u2009\u222a\u2009j+(t)$ in Eq. (56). These two terms can be computed using the expressions given from Eq. (40) through Eq. (52). Note that this paper only focuses on reliability analysis of systems with bivariate responses. Extension of the proposed method to systems with higher number of responses will be investigated in future.

Until now we have all the equations needed for the time-dependent system reliability analysis for bivariate responses.

## Numerical Procedure

We provide a flowchart for the proposed method in Fig. 3; we also summarize the main steps as below:

*Step 1*: Initialization of parametersTransform the non-Gaussian random variables and stochastic processes into standard Gaussian random variables and stochastic processes.

*Step 2*: FORMPerform the MPP search at time instant

*t*using Eq. (18); obtain the associated reliability indices, and the derivative of reliability indices, such as $\alpha i(t),\u2009\alpha \u02d9i(t),\u2009\beta i(t),\u2009\beta \u02d9i(t),\alpha j(t),\u2009\alpha \u02d9j(t),\u2009\beta j(t),\u2009\beta \u02d9j(t)$, by applying Eqs. (20), (21), (A3), and (A6).*Step 3*: Initial reliabilityCalculate the initial component reliability using Eq. (25) and initial system reliability using Eqs. (57) and (58).

*Step 4*: Upcrossing rates and outcrossing rateCompute the component upcrossing rates using Eq. (A2) and $vij+\u2212(t)$ and $vij\u2212+(t)$ using Eq. (38); then obtain the joint upcrossing rate $vi\u2009\u222a\u2009j+(t)$.

*Step 5*: IntegrationIntegrate the upcrossing rates over $[t0,\u2009ts]$ using Eqs. (14) and (25).

*Step 6*: System reliabilityObtain the system probability of failure $pf,\u2009s(t0,\u2009ts)$ using Eq. (59) or (60).

## Numerical Examples

In this section, we use three examples to demonstrate the proposed method. They are a Daniels system [48], a function generator mechanism system, and a hydrokinetic turbine system. Each example represents an important area of applications. Example 1 is a parallel system and a structural analysis problem where both random variables and stochastic processes are involved. Example 2 is a series system and a mechanism analysis problem. Even though there are no stochastic processes in the limit-state functions in Example 2, the responses of the mechanism system are still stochastic process, because the limit-state functions are functions of time and random variables. Example 3 is a hydrokinetic turbine system where the turbine blades are subjected to nonstationary stochastic river flow load.

### Example 1: A Daniels System.

Figure 4 shows a structural system under stochastic loading.

The system consists of two bars. Due to different manufacturing precisions, the two bars have different standard deviations in their dimensions. As the two bars are exposed to corrosions, their widths and heights decrease at the rates of *k*_{1} and *k*_{2}, respectively. Each of the two bars resists a load of $P(t)/2$ until both of the two bars yield. The task is to determine the time-dependent system probabilities of failure over different time intervals up to [0, 20] years.

and $X=[a1,\u2009b1,\u2009a2,\u2009b2,\u2009\sigma b1,\u2009\sigma b2]$, and $Y(t)=[P(t)]$; $\sigma b1$ and $\sigma b2$ are the yield strengths of bars 1 and 2, respectively. The parameters in Eqs. (62) and (63) are presented in Table 1.

where $\zeta =2\u2009yr$ is the correlation length. The longer is the time interval $t2\u2212t1$, the weaker is the autocorrelation.

To evaluate the accuracy of the new method, we also performed MCS using a large sample size of 10^{7}. We compared the upcrossing rates $v12+\u2212(t)$, $v12\u2212+(t)$, and $v1\u2009\u222a\u20092+(t)$ obtained from the proposed method and MCS as well.

Figures 5–7 depict the upcrossing rates $v12+\u2212(t)$, $v12\u2212+(t)$, and outcrossing rate $v1\u2009\u222a\u20092+(t)$ from both the new method and MCS.

Note that the curves of upcrossing rates and outcrossing rate from MCS are not smooth. The noise comes from the numerical discretization of stochastic process. Nevertheless, the results show the good consistency between the MCS results and those from the proposed method. This example indicates that the proposed method can produce accurate joint upcrossing rates and outcrossing rate that are needed for time-dependent system reliability analysis.

Using the outcrossing rate $v1\u2009\u222a\u20092+(t)$, we obtained the system reliability analysis result. The joint probability $pf,1\u222a2(t0,\u2009ts)$ and time-dependent system probability of failure $pf,\u2009s(t0,\u2009ts)$ are depicted in Figs. 8 and 9, respectively. The latter is also given in Table 2.

As shown in Fig. 9, the error of the new method becomes larger with a longer period of time or with a larger probability of failure. The error resource is mainly the assumption of independent crossings. It is the intrinsic drawback of the upcrossing and outcrossing rate method [33].

To show the computational cost of the new method, we also provide the numbers of function calls in Table 3. The results indicate that the new method is much more efficient than MCS.

### Example 2: A Function Generator Mechanism System.

A function generator mechanism is a mechanism used to realize a desired motion [20,21]. Such a system is shown in Fig. 10. This system consists of two function generator mechanisms. Mechanism 1, a four-bar linkage mechanism with links $B1$, $B2$, $B3$, and $B4$, generates a sine function, while mechanism 2, the other four-bar linkage mechanism with links $B1$, $B5$, $B6$, and $B7$, generates a logarithm function.

where $B=[B1,B2,\u2026,B7]$.

where $D\eta =2B7(B1\u2212B5\u2009cos\u2009\theta )$, $E\eta =\u22122B5B7\u2009sin\u2009\theta $, and $F\eta =B12+B52+B72\u2212B62\u22122B1B5\u2009cos\u2009\theta $.

In this problem, the time factor is the input angle $\theta $. There are no stochastic processes in the input variables. The vector of random variables is therefore $X=B=[B1,B2,\u2026,B7]$, and the vector of stochastic processes $Y$ is empty. Since the time factor $\theta $ appears in both functions of the motion errors, the motion errors are still stochastic processes. The motion errors should not be large, and their allowable values are denoted by $e1$ and $e2$. All the parameters are given in Table 4.

Figures 11–13 show the results of joint upcrossing rates $v12+\u2212(\theta )$, $v12\u2212+(\theta )$, and outcrossing rate $v1\u222a2+(\theta )$, from the proposed method and MCS. The sample size of MCS is 10^{7}.

The results show that the proposed method is able to estimate the joint upcrossing rate with good accuracy. Based on the joint upcrossing rates, we obtained the time-dependent system probability of failure as presented in Fig. 14 and Table 5.

The results show that the accuracy of the proposed method is good. Table 6 gives the number of function calls required by the new method and MCS.

### Example 3: A Hydrokinetic Turbine System

#### Problem Statement.

A hydrokinetic turbine system is employed as our third example [24,49]. This system is used to extract energy from river water flow, and it is subjected to stochastic flow loads during its operation. A three-dimensional model of the hydrokinetic turbine blade is created as shown in Fig. 15. The blade is 1-m long and made of steel. The turbine blade is twisted and has variable chord length along the radial direction. The hydrofoil of the blade is NREL S809. The lift and drag coefficients of the hydrofoil are available in Ref. [50].

*i*th blade, respectively, and $Y(t)=[Vr(t)]$ is the stochastic river flow velocity. Table 7 gives the statistical information of the variables involved. The mean $\mu V(t)$, standard deviation $\sigma V(t)$, and correlation function $\rho V(t1,\u2009t2)$ of the stochastic river velocity $Vr(t)$ are given by

The above information of the river flow velocity is obtained by analyzing the historical data of the Missouri river [51].

Next, we first discuss how to perform FEA analysis to obtain the maximum strain for a given river flow velocity and realization of random variables. We then present the process and result of the time-dependent reliability analysis.

#### Finite Element Analysis.

From the load analysis, it is found that the turbine blade is subjected to an edgewise moment generated from the edgewise force $FT$ and a flapwise moment generated from the flapwise force $FN$. In order to compute $FT$ and $FN$, the turbine blade is divided into 48 stations along the radial direction. After the discretization, the edgewise force $FT,\u2009i$ and flapwise force $FN,\u2009i$ at station *i* are computed using the blade element momentum theory [52] based on the geometry of the turbine blade, local pitch, and the river velocity at the station. More details about the load analysis of the turbine blade is available in Ref. [52]. After the forces at stations of the turbine blade are obtained, they are input into FEA to get the stress response of the blade. Figure 16 gives the flowchart of the stress analysis. Figure 17 plots a snapshot of the stress analysis results of the blade. From the stress response, we obtain the maximum stress and also the maximum strain.

#### Time-Dependent Reliability Analysis.

Figures 18–20 show the results of $v12+\u2212(\theta )$, $v12\u2212+(\theta )$, and $v1\u222a2+(\theta )$, from the proposed method and MCS. The sample size of MCS is 10^{6}.

The results show that the proposed method is able to estimate the joint upcrossing rate accurately even for problems with nonstationary Gaussian process. Based on the joint upcrossing rates, we obtain the time-dependent system probability of failure as presented in Fig. 21 and Table 8. The results show that the accuracy of the proposed method is good. Table 9 gives the number of function calls required by the new method and MCS.

## Conclusion

Time-dependent system reliability analysis plays a vital role in the system level optimization, lifecycle cost estimation, and decision making on maintenance and warranty. With the availability of computational models, predicting how the system reliability changes with time is possible. Making such a prediction both accurate and efficient is critical.

In this work, we proposed a time-dependent reliability method for a system with two response variables that are functions of random variables and stochastic processes. The method is based on the FORM and the upcrossing rate method (the Rice's formula). The new method can be applied to general problems with random variables, stochastic processes, and time, because it can be extended to systems with more than two response variables. With the use of FORM, the proposed method is also efficient. However, the error of FORM may affect the accuracy of the proposed method. The proposed method is therefore limited to problems when FORM is applicable.

As an upcrossing rate method, the new method may produce a larger error if the probability of failure is larger. The reason is that when the probability of failure is large, the dependency between upcrossings may become strong.

The future research based on this work may be: (1) the improvement of the accuracy, (2) the extension of the method to systems with more than two response variables, and (3) the integration of the method with optimization so that the time-dependent system reliability-based design can be performed.

### Appendix A: Time-Dependent Component Probability of Failure, $pf,\u2009k(t0,\u2009ts)$

where $\Phi (\u22c5)$ is the CDF of a standard normal random variable.

and $\varphi (\u22c5)$ is the PDF of a standard normal random variable.

*l*.

### Appendix B: Transformation of $vij+\u2212(t)$

where $fLj|Li=\beta i(t)(\u22c5)$ is the conditional PDF of $Lj$ given $Li=\beta i(t)$, and $fL\u02d9i|Li=\beta i(t),\u2009Lj=lj(\u22c5)$ is the conditional PDF of $L\u02d9i$ given $Li=\beta i(t)$ and $Lj=lj$.

where $Z\u02d9L=(\beta \u02d9i(t)\u2212\mu L\u02d9i|Li=\beta i(t),\u2009Lj=lj)/\sigma L\u02d9i|Li=\beta i(t),\u2009Lj=lj$.

### Appendix C: Derivation of $vij++(t)$

*t*, we have $Zi(t+\Delta t)\u2248zi(t)+z\u02d9i(t)\Delta t>0$ and $Zj(t+\Delta t)\u2248zj(t)+z\u02d9j(t)\Delta t>0$, then Eq. (55) becomes

## Acknowledgment

This material is based upon the work supported in part by the Intelligent Systems Center at the Missouri University of Science and Technology.

## Funding Data

Division of Civil, Mechanical and Manufacturing Innovation, National Science Foundation (Grant Nos. CMMI 1234855 and CMMI 1727329).