Time-variant reliability analysis of a continuous system with strength deterioration based on subset simulation

To conduct a reliability analysis for mechanical components, it is necessary to consider the combined influence of strength deterioration and dynamic loads. An efficient method based on subset simulation is proposed in this paper to analyze time-variant reliability by considering the strength deterioration of mechanical components in a continuous system. A gamma process is used to describe the deterioration of system strength. A model for time-variant reliability considering strength deterioration is constructed for a continuous system. A representative example and tubular cantilever structure are assessed to demonstrate the efficiency and accuracy of the proposed method. The reliability probability examples were analyzed using a first-order reliability method and benchmark results for the proposed method were derived using direct Monte Carlo simulation (MCS). The results of the proposed method and MCS are consistent, indicating that the proposed method is an effective reliability analysis method for evaluating small failure probabilities in a continuous system subjected to strength deterioration and dynamic loads.


Introduction
The reliability analysis of mechanical components is critical for improving the operational reliability and efficiency of mechanical equipment. Considering the influence of changes in the environment and uneven material surfaces, the structural and performance parameters of mechanical components are often uncertain. Some parameters are timeinvariant random variables subjected to a particular distribution (e.g., normal distribution), while others are timevariant random variables that change regularly over time or follow a random process. According to different typologies of uncertainty, analysis methods for reliability can be divided into two main categories: time-variant and timeinvariant methods. Over the past few years, several studies on reliability methods have been reported. Such methods have been widely utilized in engineering applications [1][2][3]. Castaldo et al. [1] proposed a sectional approach based on Monte Carlo simulation (MCS) to estimate the expected lifetime of a deteriorating reinforced concrete bridge pier. Huang et al. [2] proposed a new second-order reliability method with saddle point approximation for accurate, convergent, and computationally efficient estimates of the probability of failure. Zhu et al. [3] incorporated finite element simulations with Latin hypercube sampling to assess the fatigue reliability of high-pressure turbine discs under complex loadings coupled with multisource uncertainties.
In most cases, mechanical components are subjected to strength deterioration and the failure process is gradual. Several methods have been reported for analyzing component deterioration. Mori and Ellingwood [4] developed methods to evaluate the time-dependent reliability of reinforced or pre-stressed concrete structures using structural reliability principles. Li [5] postulated that the timedependent reliability analysis of a deteriorating structural system was based on the reliability of a given failure sequence for the structure. Ciampoli [6] formulated a probabilistic method based on a stochastic differential equation for the reliability assessment of structural components subjected to deterioration. Li et al. [7] developed an improved method for evaluating the time-dependent reliability of structures to assess the effects of nonstationarity in the load and resistance deterioration processes on safety and serviceability quantitatively.
Mechanical components typically bear dynamic loads during regular system operation. Therefore, it is essential to evaluate the continuous reliability of mechanical components by combining dynamic loads with strength deterioration. Rice [8] presented a novel first-passage probability formula, which has been widely used to evaluate timevariant reliability. Andrieu-Renaud et al. [9] developed a method known as PHI2 based on an out-crossing approach that solved reliability problems using classical time-invariant reliability tools, such as the first-order reliability method (FORM) and second-order reliability method. Zhang et al. [10] proposed an approach to translate timevariant reliability into static reliability by discretizing the stochastic processes of resistance and stress at calculation time. Li and Kiureghian [11] developed a method based on the principles of optimal linear estimation theory for the efficient discretization of random fields to simulate stochastic processes.
Subset simulation (SS) was originally proposed by Au and Beck [12] to estimate the small failure probability of high-dimensional parameter space problems. Small failure probability is expressed as the product of the large conditional probabilities of selected intermediate failure events. This method has been widely used in statistics, applied mathematics, and structural engineering. Au and Beck [13] used SS to effectively simulate procedure for seismic performance assessment of structures in the context of modern performance-based design. Vahdatirad et al. [14] estimated low-event probabilities using SS to analyze the first natural frequency of a turbine supported by a surface footing. Norouzi and Nikolaidis [15] re-evaluated the reliability of a dynamic system subjected to stationary Gaussian stochastic excitation for various load spectra by reweighting the results of a single SS. Song et al. [16] analyzed the reliability sensitivity of failure probabilities by transforming the distribution parameters of basic variables into a set of conditional failure probabilities using SS. Bourinet et al. [17] presented an approach referred as 2SMART for estimating small failure probabilities by considering SS from the perspective of support vector machine classification. Zuev et al. [18] presented a modified Metropolis algorithm based on Markov chain MCS and used SS for sampling from conditional distributions to compute small failure probabilities in general high-dimensional reliability problems.
Recently, SS has been widely used to evaluate dynamic structural reliability [19][20][21][22]. Li et al. [19] developed a generalized SS approach to estimate the failure probabilities of multiple stochastic responses. Wang et al. [20] presented an improved SS method with a splitting approach to estimate the time-dependent reliability of systems with random parameters excited by stochastic processes. Yu and Wang [21] proposed a novel time-variant reliability analysis method using failure process decomposition to transform time-variant reliability problems into time-invariant problems for dynamic structures with uncertainty. Yu et al. [22] proposed a novel time-variant reliability analysis method for multiple failure modes and temporal parameters based on a combination of the extreme value moment method and improved maximum entropy method. However, time-variant reliability problems with strength deterioration for mechanical components have rarely been evaluated using SS. In this paper, a novel approach is proposed for the continuous system reliability analysis of mechanical components using SS. Expressions for the performance functions of a continuous system are constructed using a model for time-variant reliability with strength deterioration and the failure probabilities of mechanical components are derived.
The remainder of this paper is organized as follows. Section 2 describes a model for deterioration processes and presents a performance function for a continuous system, as well as a model of time-variant reliability with strength deterioration. Section 3 describes an approach for evaluating continuous system reliability using SS. Section 4 presents two numerical examples to demonstrate the efficiency and accuracy of the proposed approach. Conclusions are provided in Section 5.

Model of time-variant reliability with strength deterioration
The gamma process proposed by Abdel-Hameed [23] is the most appropriate process for describing the deterioration of strength. Excluding cyclic loads, a given system is largely influenced by random environmental factors, such as random vibrations in the external environment, stochastic variation in ambient temperature, or corrosion. These factors lead to a monotonic decrease in strength and stochastic process evolution over time. Furthermore, the gamma process is a continuous-time stochastic process with Time-variant reliability analysis of a continuous system with strength deterioration based on … 189 independent, nonnegative increments (e.g., variation in tubular cantilever structural strength). Therefore, it is suitable for modeling gradual damage that accumulates monotonically over time in a sequence of small increments, such as wear, fatigue, corrosion, creep, and a degrading health index. Additionally, relatively straightforward mathematical calculations are a major advantage of modeling deterioration processes as gamma processes. Van Noortwijk [24] described the statistical properties of gamma processes as a deterioration model. The symbol Y(t) is used to denote the deterioration at a time t, where t ! 0. Therefore, according to the definition of a nonstationary gamma process, the probability density function (PDF) of Y(t) can be expressed as f YðtÞ ðyÞ ¼ GðyjtðtÞ; uÞ where I (0,?) (y) = 1 for y [ (0,?), I (0,?) (y) = 0 for y (0,?), and CðtÞ ¼ R 1 0 y tÀ1 e Ày dy is the gamma function for a positive shape parameter t. Furthermore, the gamma process is a stochastic process with a shape function t(t), which is non-decreasing and right-continuous, and positive scale parameter u. Therefore, the expected value can be expressed as The variance can be written as varðYðtÞÞ ¼ Experimental studies have confirmed that the expectation of deterioration at a time t is generally proportional to the following power law where a, b, and c are positive physical constants. For applying the gamma process model to engineering problems, statistical methods are required to estimate the parameters of the shape function t(t) = ct b and scale parameter u. Through various calculations, practical knowledge [25][26][27][28] can be obtained regarding the shape function in terms of the parameter b in Eq.(4). Consequently, the value of b is assumed to be known, but the values of c and u are unknown. The estimators of c and u can be obtained using the most common method of parameter estimation, namely the maximum likelihood. Cinlar et al. [26] provided a method for transforming a nonstationary gamma process into a stationary gamma process. The corresponding observations of cumulative deterioration fy i t j i t ¼ 1; 2; Á Á Á ; n t g, where y 0 \y 1 \y 2 \ Á Á Á \y n t and y 0 = 0, are obtained at inspection time points ft i t ji t ¼ 1; 2; Á Á Á ; n t g, where t 0 \t 1 \t 2 \ Á Á Á \t n t and t 0 = 0. The likelihood function for the observed deterioration increments Based on the above analysis, the residual strength S(t), when formulated as a stochastic process, is equal to the initial strength S 0 , and is assumed to decrease monotonically based on a deterioration function Y(t). Therefore, the general form of assumed residual strength can be expressed as where S 0 is a random variable with a normal distribution, which is closely related to component elements, internal defects, heat treatment, and work-hardening of materials. Y(t) is a stochastic process depending on the conditions of suffered stress and operating environment. Therefore, S 0 and Y(t) are uncorrelated. Based on the above definitions and the stress-strength interference theory, the performance function for timevariant reliability with strength deterioration for a continuous system can be expressed as where the stochastic process X t ð Þ is composed of the initial strength S 0 , deterioration function Y(t), and external loads Z t ð Þ which include the static and dynamic loads at the time t. Therefore, max rðZ; tÞ ð Þdenotes the maximum dynamic stress magnitude imposed by external loading during the time interval.
Consequently, the failure probability of a continuous system with strength deterioration can be calculated as The model of time-variant reliability can be expressed as MCS is generally a feasible method to estimate the probability of failure once a system performance function is formulated. However, MCS is not appropriate for obtaining small probabilities of failure (e.g., P F \10 À3 ) because it suffers from a crucial lack of efficiency to achieve the required accuracy. Therefore, in this study, a Markov chain Monte Carlo (MCMC) method based on the Metropolis algorithm was used to efficiently compute small failure probabilities in a continuous system bearing dynamic loads and suffering from gradual failure caused by strength deterioration.

SS-based time-variant reliability analysis with strength deterioration
The core idea of the SS method proposed by Au and Beck [12] is to estimate the frequency of a rare event based on the frequencies of more common events in a sequence of intermediate failure events. SS has been widely used to analyze reliability effectively. However, solving the problem of time-variant reliability analysis with strength deterioration using SS has received little attention. Based on the performance function described in Section 2, SS is expanded in this section to handle the problem of timevariant reliability analysis with strength deterioration for a continuous system.

SS for a continuous system
Regardless of whether or not the stochastic process X t Þis defined with respect to continuous-time or discrete-time variables, only the values for a process at a discrete number of time instances are generated. For the sake of convenient calculation, this time grid is defined by 0; t 1 ; t 2 ; Á Á Á ; t n t À1 ; t n t , where t i t ¼ ði t =n t Þt for t [ 0 and i t ¼ 0; 1; Á Á Á ; n. Then, X t ð Þ can be expressed as vectors of random variables X 1 ; X 2 ; Á Á Á ; X n t À1 ; X n t (i.e., X i t j i t ¼ 1; f 2; Á Á Á ; n t g).
The symbols fF i t ji t ¼ 1; 2; Á Á Á ; n t g define the continuous system failure events at the time instances t i t . Furthermore, Consequently, the probability of failure at a time t i t can be expressed as where where n is the number of variables. The symbol q j ðÁÞ represents the univariate PDF for each component of X i t . Equation (9) calculates the product of a sequence of intermediate conditional probabilities P F i t ;i F i t ;iÀ1 À Á i ¼ È 2; 3; Á Á Á ; mg and the first level PðF i t ;1 Þ. To make the conditional probabilities in Eq.(9) sufficiently large to be estimated efficiently, it is important to select intermediate failure events appropriately. Particularly, the prior determination of b i is not a trivial task. Therefore, in this study, b i was selected such that the estimated conditional probabilities were equal to a fixed value p 0 2 ð0; 1Þ.
This was accomplished by setting the intermediate threshold values sorted in an ascending order. The symbols Additionally, the probabilities PðF i t ;1 Þ and É are required to obtain P F;t i t from Eq. (9). MCS can be used to estimate PðF i t ;1 Þ as where Similarly, P F i t ;i F i t ;iÀ1 À Á can be computed using a formula similar to Eq.(10). However, to obtain samples according to the conditional distributions q efficiently, simulation should be carried out using Markov chain MCS based on the modified Metropolis algorithm proposed by Papaioannou et al. [29], rather than traditional MCS.

Implementation procedure
In general, SS for a continuous system involves six main steps, which are detailed below.
(i) Discretize the duration of continuous system operation T using an appropriate sampling interval Dt to obtain a number of time instances n t ¼ T=Dt and discrete time instances (ii) At time instance t 1 , generate N l i.i.d. samples X P g X 1 ð Þ\b 1 f gequal to p 0 , assign the (p 0 N l ? 1)th value to b 1 . (iv) For every i ¼ 2; 3; Á Á Á, generate 1 À p 0 ð ÞN l extra conditional samples based on the p 0 N l seed samples in the failure region F 1;iÀ1 ¼ g X 1 ð Þ\b iÀ1 f g using the modified Metropolis algorithm [29,30] to obtain a total of N l samples (i.e., X 1;k l ) in F 1;iÀ1 ¼ fgðX 1 Þ\b iÀ1 g. Similarly, calculate the values of the performance function using these N l samples and sort gðX ðiÀ1Þ 1;k l Þ k l ¼ 1; 2; Á Á Á ; N l n o in ascending order again. Assuming that g are equal to p 0 , take the (p 0 N l ? 1)th values of the sequence as b i to obtain 3; Á Á Á ; m À 1g, and Eq.(11) into Eq.(9), the failure probability at time instance t 1 can be denoted as (vi) Repeat steps (ii)-(v) for the time instances t i t ; i t ¼ 2; 3; Á Á Á ; n t . Then, the continuous system failure probability can be denoted as In this study, SS was used to calculate the time-variant reliability for a continuous system with strength deterioration. The small failure probability is transformed into the product of larger probabilities, thereby achieving a reduction in the number of required calculations and improving calculation efficiency while maintaining high precision. Figure 1 presents the implementation process of the proposed method concisely.

Example applications
In this section, two representative examples of a high-dimensional stochastic problem and continuous engineering system are used to demonstrate the accuracy and efficiency of the proposed method. MCS is the most accurate method for solving the problem of time-variant reliability analysis. Considering the nonlinearity of the performance function, the FORM [31] was adopted based on its excellent performance and simplicity. The results obtained using the proposed method were compared with those obtained using the FORM and MCS. To achieve the desired accuracy, MCS requires 10 kþ2 À 10 kþ3 samples if the failure probability is on the order of 10 -k [32]. The variance values for the random parameters can be obtained through reliability testing or statistical analysis of experimental data, or can be estimated based on the following principles if there is no relevant data. The variance can be estimated based on a variation coefficient. In general, for mechanical performance parameters related to metal materials, this coefficient is set to 0.05 [33]. If the system is subjected to the influence of a large number of independent factors (no dominant factor), random variables commonly follow a normal distribution [34].

High-dimensional example
This example was chosen to assess the performance of the proposed method for handling relatively high numbers of random variables. It is known that increasing the dimensionality of the random variable space introduces various challenges (e.g., curse of dimensionality [35]). We propose analyzing the academic example of one failure mode of a bevel gear, where z 12 is the diameter of the small gear and z 13 is the length of the contact line at the midpoint of the tooth surface. z i i ¼ 1; 2; Á Á Á ; 11; 14; Á Á Á ; 17 ð Þ is the coefficient of the bevel gear and F is the nominal tangential force on a reference circle at the midpoint of the tooth width. The statistical distributions and parameters for the random variables are listed in Table 1. According to international standards [36], the contact stress of a bevel gear rðZ; tÞ is expressed by the load on the bevel gear transmission. The deterioration process for the contact fatigue strength of the tooth surface S(t) can be described by a gamma process. The performance function is written as where z i ji ¼ 1; 2; Á Á Á ; 17 f g , and F and S 0 are all independent. Here, u = 1.611 1. Based on the S-N curves, the parameters for the deterioration process of system strength are u = 1.382 2 9 10 -8 , b = 0.05, and c = 4.975 1 9 10 10 .
According to the performance function, the failure probabilities of the continuous system P F were evaluated using the proposed method, FORM, and MCS. The results are listed in Table 2. Figure 2 presents the cumulative probability distributions of the continuous system as calculated by the proposed method and MCS at each time instance on logarithmic scales. One can see a clear connection between the cumulative probability, value of the performance function, and time of service. As shown in Table 2 and Fig. 2, the failure probabilities and cumulative probabilities of the continuous system estimated by the proposed method agree well with those acquired from MCS. Compared to the results of the proposed method, the failure probabilities for the continuous system obtained by the FORM are significantly higher or lower than those evaluated by MCS.
To illustrate the high efficiency of the proposed method relative to MCS, Figs. 3 and 4 compare the coefficient of variation (COV) of the failure probabilities d i t [12,37] and

Tubular cantilever structure
The tubular cantilever structure presented in Fig. 5 was considered as a second example. This example was studied by Madsen et al. [38], and Du and Chen [39] as an engineering design problem. The cantilever structure was analyzed based on its resistance to yielding caused by bending and sheering stress. Seven variables were involved in this example: two stochastic process variables F 1 (t) and T(t), and five random variables S 0 , F 2 , F 3 , d, and h. S 0 was the initial strength of the tubular cantilever structure and strength deterioration was modeled as a gamma process. F 1 (t), F 2 , F 3 , and T(t) are the three external forces and torque imposed on the cantilever structure, respectively, and d and h are the dimensions of the cross section. Their statistical distributions and parameters are listed in Table 3 and are all independent. Additionally, the locations of the force points and angles relative to vertical of the forces  F 1 (t) and F 2 were L 1 = 60 mm, L 2 = 120 mm, h 1 = 5°, and h 2 = 10°. As shown in Fig. 5, F 1 (t), F 2 , F 3 , and T(t) cause the cantilever structure to experience bending and torsion.
Based on analysis of mechanical model, the bending and sheering stresses can be calculated as follows s xz ðtÞ ¼ TðtÞd According to Eqs. (14) and (15), the maximum stress can be expressed as Based on the deterioration of strength, the performance function takes on the following form gðX; tÞ ¼ SðtÞ À max rðZ; tÞ ð Þ ¼ S 0 À YðtÞ À r max ðtÞ; ð17Þ where Y(t) is defined by parameters u = 1.486 3 9 10 -6 , b = 0.2, and c = 2.834 9 9 10 7 . The failure probability results are listed in Table 4. These results were obtained by the proposed method, FORM, and MCS based on Eq. (17). A comparison of the cumulative probability distributions of the tubular cantilever structure failure system as calculated by the proposed method and MCS at each service time instance is presented in Fig. 6. Figures 7 and 8 indicate that the computational time required by the proposed method is still much shorter than that required by MCS in the case of For the instance of t = 1, the number of samples required for MCS was 10 6 , but the proposed method required only 10 3 initial samples. The computational time required by the proposed method (t s = 6.951 9 s) is    significantly shorter than that required for MCS (t s = 338.288 5 s). And the P F,1 value estimated by the proposed method is 3.81 9 10 -4 , which agrees well with the result obtained from MCS, 3.80 9 10 -4 . However, the P F,1 value calculated by the FORM is 9.87 9 10 -6 , which differs significantly from the MCS result. This discrepancy is largely based on the errors induced by the nonlinearity of the performance function for the tubular cantilever structure system and variety of the distributions of random variables. Additionally, the failure probabilities of the events were all small. These results further demonstrate that SS modified for a continuous system is an accurate and

Conclusions
A novel reliability analysis method was proposed in this paper. The proposed method combines the advantages of SS with the discretization of continuous systems. Continuous service times are divided into appropriate time intervals to determine the variation in failure probabilities for continuous systems at each time instance using SS. This method evaluates the small failure probabilities of timevariant reliability with strength deterioration efficiently and accurately. The proposed model for time-variant reliability with strength deterioration performs well on reliability analysis problems for continuous systems subjected to strength deterioration and dynamic loads. Two examples demonstrated that the proposed method estimates P F;i t values reasonably accurately compared to common reliability analysis methods (e.g., the FORM), particularly if the failure probabilities are as low as 10 -4 . Furthermore, the computational time required for the proposed method is considerably shorter than that required for MCS. The proposed method can be utilized for continuous systems with small failure probabilities as an effective reliability analysis method. Additionally, it provides a new perspective for time-variant reliability problems with strength deterioration characterized by small failure probabilities, multiple random variables, and nonlinearity of performance functions.