Time-dependent reliability analysis of corroded RC beams based on the phase-type fitting method

Highlights Abstract ▪ A novel phase-type (PH) fitting method is developed for time-dependent reliability. ▪ A reliability model of reinforced concrete beams is formulated by considering the time-dependent chloride diffusion coefficient.


Introduction
The corrosion of steel bars in concrete stands as a prominent determinant leading to the failure of reinforced concrete (RC) structures.The resultant damage to these structures can be categorized into two distinct types: concrete carbonation and chloride penetration [3,18,42,53].However, in practical engineering, the rate of chloride ion penetration surpasses that of concrete carbonation [63].Hence, the chloride penetration has stronger impact on reinforcement corrosion and structural deterioration.
Two primary sources of chloride ions have been identified [22]: the first source involves the addition of chloride during the mixing and pouring phase, including calcium chloride and sodium chloride.The second source originates from the external environment during the service period, encompassing deicing salts used in winter, seawater, and sea breeze.Chlorides infiltrate the concrete through diffusion, driven by concentration gradients.Once the chloride concentration at the depth of a steel bar exceeds the chloride threshold level, the protective passivation film on the surface of the steel bar deteriorates due to the combined effects of water, oxygen, and chloride [40,46,50].Subsequently, corrosion occurs within the bar.The corrosion of steel reinforcements can lead to concrete fractures through cracking, delamination, and spalling of the concrete cover.Moreover, it results in a reduction in the crosssectional area of both the concrete and reinforcements, significantly compromising the serviceability, strength, safety, and lifespan of concrete structures [44].
Fick's second law provides an effective description of chloride penetration in concrete.In 1972, Collepardi et al. [10] introduced a method for calculating the chloride diffusion coefficient based on Fick's second law, which subsequently gained widespread adoption [4,8,9,13,36,56].The chloride diffusion coefficient plays a pivotal role in the transportation of chlorides within concrete.Many diffusion models, rooted in Fick's second law, assume a constant value for the chloride diffusion coefficient [2,10,18,44,23].However, in practical scenarios, particularly in marine environments, the diffusion of chloride ions into concrete exhibits nonlinear and timedependent behavior [24,61].Consequently, the actual chloride diffusion coefficient varies over time, rather than remaining constant.
The performance of RC structures is influenced by various factors, including material properties and environmental conditions.Given the material heterogeneity, environmental complexity, and uncertainties associated with both, a probabilistic approach is necessary to study the corrosion damage of RC structures [62].Additionally, when considering the time-dependent chloride diffusion coefficient and stochastic characteristics of loads, it is essential to perform timedependent reliability analysis of RC structures.Saassouh and Lounis [44] used the first-order reliability method (FORM) and second-order reliability method (SORM) to predict the timedependent probability of a RC highway bridge deck.El Hajj Chehade et al. [14] proposed a simplified time-dependent reliability model for a set of 21 simply supported RC T-beam bridges under variable traffic load considering the concrete creep and shrinkage.Some other related studies include Bagheri et al. [5], Guo et al. [20], Fan et al. [17] and Yang and Li [65].
Reliability is quantified as the integration of the irregular region in multidimensional space, and direct integration methods exhibit limited feasibility.For intricate reliability analysis problems, Monte Carlo simulation (MCS) [7,67] stands as a potent and versatile technique.Its advantage lies in obtaining precise numerical solutions directly based on simulation test results.However, MCS entails a significant computational burden, which poses a major drawback.The FORM [66] and SORM [71] offer a relatively lower computational cost compared to MCS for evaluating reliability.
Nonetheless, as the nonlinearity of the limit state function (LSF) increases, the approximate accuracy of these methods cannot be guaranteed [29].
The phase-type (PH) fitting method presents a substantial reduction in computational cost while ensuring reliability accuracy.Furthermore, this method possesses the advantageous property that any positive dataset can be approximated as a PH distribution using the Expectation Maximization (EM) algorithm.As a result, the PH fitting method has gained widespread usage in the field of reliability engineering [1,12,33,57,58].However, previous studies have relied on trial methods [31] to determine relevant parameters, which involve a complex, inefficient, and unscientific operational process.
Different from the original PH fitting method [31,43,55], the proposed novel method identifies the parameters in an automatic, efficient, and scientific iterative manner.The stopping criterion is adopted to simplify the process of parameter estimation.
Considering the limitations of the aforementioned research, this paper focuses on the time-dependent reliability analysis of corroded RC beams using the PH fitting method.The paper aims to achieve three primary objectives: 1) developing a reliability model for RC structures that incorporates the timedependent chloride diffusion coefficient; 2) proposing a novel approach to obtain the parameters of the PH fitting method in a simple, efficient, and scientifically rigorous manner; 3) demonstrating the accurate and efficient application of the novel PH fitting method in handling time-dependent reliability analysis of corroded RC beams.
The paper is organized as follows.Section 2 establishes the framework for the time-dependent reliability analysis of corroded RC beams.Section 3 introduces the novel PH fitting method.In Section 4, the proposed methodology is validated through numerical examples and a case study involving corroded RC beams.Finally, Section 5 presents the conclusions derived from this study.

Review of time-dependent reliability analysis
In engineering applications, a multitude of uncertainties arises throughout the life cycle of structures, spanning the stages of design, construction, and operation.These uncertainties are commonly represented by random variables or stochastic processes.In the context of time-dependent reliability analysis, the general form of a LSF can be expressed as G(X, Y(t), t), where X represents a vector of random variables, Y(t) denotes a vector of stochastic processes, and t signifies the time parameter.Over a specified time interval [0, TL], failure occurs if there exists a specific time instant t within [0, TL] where the LSF value falls below zero [26]: As a result, the time-dependent probability of failure Pf (0, T) over a time period [0, T] (0≤T≤TL) is a monotonically increasing function with respect to T. It can be defined as: And time-dependent reliability is given by: (0, ) = (∀ ∈ [0, ], (, (), ) > 0) , 0 ≤  ≤   = 1 −   (0, ), 0 ≤  ≤   (3)

Modeling of time-dependent reliability analysis for corroded RC beams
To simplify the problem, only the corrosion model for steel reinforcement bars is considered, without considering the corrosion model for concrete.This research is the foundation of our further works when the corrosion model for concrete is also taken into account.
The primary focus of this paper revolves around the corrosion phenomenon in RC beams caused by chloride ion attack.To describe the process of chloride penetration in concrete, Fick's second law is employed [44]: where x is the depth from the concrete surface; C(x, t) is the concentration of chlorides at depth x and time t; D is the chloride diffusion coefficient.
In the case of a constant chloride diffusion coefficient D, and considering the following initial and boundary conditions, the Crank's solution can be derived [11]: )] ( 6) where Cs is the surface chloride concentration, erf () is the error function or Gauss error function.
In practical scenarios, the chloride diffusion coefficient may vary with time.For modeling the time-dependent chloride diffusion coefficient, a power law relationship is commonly recommended [32]: With the cement hydration proceeding, the concrete pore structure is refined and the connectivity of pores significantly decrease.As the result, the chloride diffusion coefficient decreases with time [48].Based on the data from long-term field and laboratory studies, a power law relationship is proposed by Thomas and Bamforth [54], and the function is conmmonly adopted in the later researches [32,59].
where m is the age factor and Dref is the chloride diffusion coefficient at reference time tref =28 days.However, the defect of this model is that the historic change of the diffusion coefficient for a given exposure duration is not considered, which may lead to an erroneous judgment on chloride ingresses and thus the risk of chloride-induced corrosion [41].
To estimate the crucial parameter m, one can refer to the experimental study conducted by Mangat and Molloy [35]: where W/C denotes water cement ratio.It appears that the value of m directly depends on the mix proportion [34], which is represented by W/C.The Eq. ( 9) is derived by a linear regression analysis with limited data of m and W/C, and the accuracy of predicted relationship can be improved using more data [35].
According to reference [44], the serviceability limit state imposes restrictions on the normal use of a structure, which may involve excessive deformation, vibration, and localized damage such as cracking, spalling, and corrosion.It is widely recognized that when the chloride concentration at the steel location exceeds the chloride threshold level, corrosion of the reinforcement occurs, subsequently leading to cracking and spalling of the concrete cover.Hence, the LSF for corrosion initiation can be formulated as [68,69]: where Ccr is the chloride threshold level; c represents the depth of the concrete cover over the steel; C(c, t) denotes the concentration of chlorides at depth c and time t.
By considering the LSF as G(X, Y(Tin), Tin)=0, the time at which corrosion initiation occurs, denoted as Tin, can be determined by solving the equation: Since parameters Ccr, Cs, c, Dref are all random variables, Tin is also a random variable, and the LSF can be rewritten as: Utilizing the MCS method, the time-dependent reliability can be computed based on the aforementioned LSF.Notably, the random variable Tin and the time parameter t are completely separated in Eq. (12).By substituting Eq. ( 12) into Eq.( 2), the time-dependent probability of failure can be determined as: It is evident that the essence of the above formula lies in the cumulative distribution function (CDF) of Tin.In other words, the time-dependent probability of failure Pf (0, T) can be evaluated by analyzing FTin(T).

A novel PH fitting method
In view of the limitations inherent in existing reliability analysis methods (such as MCS, FORM, SORM, etc.), it is necessary to address the challenges of accuracy and efficiency.To overcome these challenges, a novel PH fitting method is proposed to accurately and efficiently handle time-dependent reliability problems.
The PH fitting procedure involves estimating the parameters of a PH distribution based on sample data or with respect to another known distribution [55].The concept of PH distributions was initially introduced by Erlang [16] and later extended by Neuts [37,38].
PH distributions can be regarded as a generalization of the exponential and Erlang distributions [45].The representations of PH distributions are depicted in Fig. 1. expressed as [64]: The steps of the novel PH fitting method are outlined as follows:

1) Initialization
Step 1: Set the initial order as N=1, and initialize the global maximum log-likelihood value as Lmax= -10 5 .
Step 3: Employ the integer splitting method to calculate the number of possible cases of state transition SN and obtain the corresponding shape parameters RN : {r1, r2, ...,   }.

3) EM algorithm
The EM algorithm is a valuable tool for estimating parameters from a given data trace [43].Thus, it is suitable to utilize the EM algorithm to determine the parameters of the HErD.The EM algorithm is executed as follows for a given .
Step 6: Update the parameters using   = ( 1 ,  2 , ⋯ ,   ,  1 ,  2 , ⋯ ,   ), expressed as Step 9: The log-likelihood values at θi (i=1, 2, ..., SN) can be calculated as And then, the current maximum log-likelihood value LN, the current optimal output parameters θ * N, the optimal index ID, and corresponding optimal r * N can be given by The flowchart of the PH fitting method is shown in Fig. 3. Fig. 3. Flowchart of the PH fitting method.

Illustrative examples
In the field of reliability, conventional distributions commonly used include the lognormal, normal, gamma, exponential, Weibull, and type I extreme value (also known as the Gumbel) distributions [6,21,25,27,28,30,39,49,51,52].In this section, these six typical distributions are employed to validate the effectiveness of the novel PH fitting method.Additionally, the proposed method is applied to an example involving corroded RC beams.

Numerical examples
In the first example, four datasets with sample sizes of K=50, 250, 500 and 10 4 are generated from the lognormal distribution, where the mean value is μ=2 and the standard deviation is σ=1.
The iterative process of the PH fitting method for this example is presented in detail in Table 1.  4 .Once the optimal parameters of the HErD are obtained, the PDF and CDF can be derived using Eqs.( 14) and ( 15).
The empirical PDF (i.e., the outer contour of the frequency histogram), the PDF obtained through kernel density estimation, as well as the PDF resulting from Gaussian fitting and PH fitting, are compared with the PDF of the lognormal distribution with μ=2 and σ=1.The results are presented in Fig. 4 a) Table 2.The optimal parameters for different distribution datasets with K=500.
Distribution dataset with K=500

An example of time-dependent reliability analysis of corroded RC beams
To further illustrate the application of the PH fitting method, an example involving corroded RC beams is presented.The LSF for this example focuses on the corrosion initiation of the reinforced steel.The parameters utilized in the example are derived from relevant literature [15,47,60,70] and are summarized in Table 3.As indicated in Eq. ( 9), the age factor m is determined by the water cement ratio (W/C).In this particular example, the value of W/C is selected as 0.35, 0.4, or 0.45.Consequently, the corresponding values of m are calculated using Eq. ( 9) as 0.275, 0.4, or 0.525, respectively.Assuming Dref =12×10 -12 m 2 /s, the relationships between the chloride diffusion coefficient D(t) and time t are plotted in Fig. 7.

Fig. 7. The variation of D(t).
It is illustrated that D(t) exhibits a rapid decrease during the early stages and a slower decline in the later period.To ensure that D(t) decreases with t [61], it is necessary to limit the value of m within the range of 0 to 1. Gjorv [19] suggested setting m to 0.4, a value consistent with the empirical formula [34].
Therefore, m=0.4 (corresponding to W/C =0.4) is adopted in this study.Subsequently, the LSF in Eq. ( 12) can be expressed as follows: (, (), ) =   −  = { In the PH fitting method, 500 random samples of W are initially generated according to their respective distributions.
Based on these random samples, 500 corrosion initiation times Tin can be computed using Eq.(25).The PDF of Tin is then approximated using kernel density estimation, Gaussian fitting, and PH fitting methods.These approximations are compared with the empirical PDF, and the results are presented in Fig. 8. Figure 8 clearly demonstrates that the corrosion initiation time Tin can be effectively approximated using the PH fitting method.The reason why the PH fitting method outperforms other methods is that the PH method can make full use of the positive dataset, accurately capture its distribution characteristics and efficiently estimate relevant parameters.
Directly approximating to some conventional distributions may ignore the features of the dataset, and result in significant errors.
The iterative process to obtain the optimal parameters for the PH fitting method is outlined in Table 4.With these optimal parameters, the CDF of Tin can be calculated using Eq.(15).According to Eq. ( 13), it is known that the timedependent probability of failure Pf (0, t) (t  [0, 100] years) can be evaluated through the CDF of Tin, i.e., FTin(t).
In the MCS [7,67] method, the time interval is discretized with Δt set as 1 year, resulting in a total of s=101 time nodes.The MCS method The PH fitting method The t-IRS method In terms of computational efficiency, it can be seen from Table 5 that the MCS method requires 101×10 5 function calls to obtain the reliability.In comparison, the proposed PH fitting method only necessitates 500 function calls, showing its efficiency in reducing the computational burden significantly.
Furthermore, when considering the reliability accuracy, the errors of the PH fitting method for all time intervals are consistently below 5%, and in some cases, even below 0.1%.
This demonstrates that the PH fitting method maintains the accuracy of reliability assessment.Compared with the MCS method, it can be concluded that the PH fitting method achieves a remarkable reduction in computational burden while ensuring the reliability accuracy.
However, the efficiency of the proposed method is lower than that of t-IRS, but its accuracy is always better than t-IRS.
The PH fitting method is a universal method for solving reliability problems with positive dataset.For more complex structural systems, the proposed method is required to incorporate with other system reliability strategies.Therefore, it is vital to accurately determine relevant statistical information of parameters and quantitatively define the parameter uncertainties.

Conclusions
This paper studies the time-dependent reliability analysis of corroded RC beams based on the PH fitting method.A reliability model of RC structures is first formulated by taking the timedependent chloride diffusion coefficient into consideration.
Next, a new strategy is cooperated with EM algorithm to simply, efficiently and scientifically obtain the parameters of PH fitting method.The effectiveness of the new method is then demonstrated using several numerical examples.Finally, the novel PH fitting method is applied to an example of corroded RC beams.The proposed method is shown to be an useful way in calculating the time-dependent reliability, and the results show excellent computational efficiency and accuracy compared with the MCS method.
The main contributions of this paper are threefold: 1) formulating a reliability model of RC structures by taking into account the time-dependent chloride diffusion coefficient; 2) a new strategy is utilized to simply, efficiently and scientifically obtain the parameters of PH fitting method; 3) using the novel PH fitting method to accurately and efficiently deal with timedependent reliability analysis of corroded RC beams.
The new method is limited to the property of PH distribution that only the positive support can be approximated as the PH distribution.In other words, the PH distribution is inapplicable to approximate the dataset with negative numbers.
Further studies will include: 1) considering the corrosion model for concrete; 2) taking into account the correlations between random variables; 3) applying the new method to more complex civil engineering problems; 4) integrating the proposed method with different optimization algorithms to optimize practical engineering problems.

Fig. 1 .
Fig. 1.Representations of PH distributions.(a) Exponential distribution.(b) Erlang distribution.(c) Hyper-Erlang distribution.As shown in Fig. 1 (a), the exponential distribution is the simplest PH distribution with an order of 1, characterized by a single parameter λ1.In Fig. 1 (b), the Erlang distribution is featured with two parameters: the rate parameter λ1 and the shape parameter r1.The state transition graph of a hyper-Erlang distribution (HErD) is presented in Fig. 1 (c).The parameters of a HErD include the number of Erlang branches M, the initial probabilities [α1, α2, …, αM], the rate parameters [λ1, λ2, …, λM], and the shape parameters [r1, r2, …, rM], where α1+α2+…+αM=1.Moreover, the order of the HErD, represented by N, is determined by the sum of the shape parameters: N=r1+r2+…+rM.When M=α1=1, the HErD degenerates into the Erlang distribution.Moreover, the exponential distribution is a special case of the HErD with M=α1=r1=1.Hence, the HErD can be considered as the generalization of the exponential and Erlang distributions.The key of the PH fitting method lies in the determining the parameters of the HErD.Once the relevant parameters are determined, the probability density function (PDF) and CDF of a hyper-Erlang random variable can be Fig. 2.

Fig. 6 .
Fig. 6.Comparison of the fitting results for different distribution datasets with K=500.(a) Normal.(b) Gamma.(c) Exponential.(d) Weibull.(e) Type I extreme value (or Gumbel).It can been seen from Fig. 5 that the PDF obtained by different methods are close enough, indicating the effectiveness of the PH fitting method.
where t represents the time parameter varying within [0, 100] years, X=[Ccr, Cs, c, Dref] T denotes the random variables.Since Y(t) does not appear in this case, one can obtain W=[X, Y(t)]=[X, Z]=[Ccr, Cs, c, Dref] T .Both the PH fitting method described in Section 3 and the MCS method are employed to calculate the time-dependent reliability.

Fig. 8 .
Fig. 8.Comparison of different approximation methods for PDF.

10 5
MCS samples of W are generated, and the corresponding values of the LSF are computed at each time node.Therefore, the MCS method requires Ncall=101×10 5 function calls in total.For the t-IRS [30] method, 30 initial samples are first generated to construct an initial instantaneous response surrogate model.And extra 54 samples are then selected to update the surrogate model.The comparison results of the three different methods are presented in Fig. 9 and Table5.

Fig. 9 .
Fig. 9. Comparison of the MCS, PH fitting and t-IRS methods.

Furthermore, the sensitivityFig. 10 .
Fig. 10.Failure probability results with different means.(a) Ccr.(b) Cs.(c) c.(d) Dref.It can be seen from the sensitivity analysis that the means of random parameters have great effects on the final reliability results.With the increase of Ccr and Dref or the decrease of Cs and c, the failure probability of the RC structure will increase.

Table 1 .
Iterative process of the PH fitting method for the lognormal distribution dataset.

Table 3 .
Distribution of random parameters for the corroded RC beams.

Table 4 .
Iterative process of corroded RC beams.

Table 5 .
Failure probability results for the MCS, PH fitting and t-IRS methods.