Time-variant reliability using the PHI2 method

Time-variant reliability using the PHI2 method : principles and validation by Monte Carlo simulation C. Andrieu-Renaud & M. Lemaire LaRAMA, IFMA/Université Blaise Pascal, 63175 Aubière, France

B. Sudret EDF R&D, MMC, Site des Renardières, 77818 Moret-sur-Loing, France

Keywords: structural reliability, time variant reliability, outcrossing approach, FORM/SORM ABSTRACT: The study of the reliability of an aging structure requires taking into account the influence of time. The classical approach for performing time-variant reliability analysis relies upon the computation of outcrossing rates (of the limit state surface) by asymptotic formulae. This method requires advanced mathematics and ad-hoc implementation. An alternative approach is possible, which allows to directly compute the outcrossing rate by means of system reliability analysis. This method is called PHI2 in this paper. Thus classical time-invariant reliability tools can be used. Both approaches are compared to closed form solution on simple examples. The quality of the computed upper bound to the probability of failure is assessed by Monte Carlo simulation. A stationary and a non-stationary application examples are proposed. The numerical features and the computational cost of each method are compared. Finally, a more realistic example is solved, namely that of the failure of a concrete beam submitted to creep and a stochastic load.

1 INTRODUCTION Structural reliability analysis aims at computing the probability of failure of a mechanical system with respect to a prescribed failure criterion by accounting for uncertainties arising in the model description (geometry, material properties) or the environment (loading). In the context of durability analysis, i.e. the assessment of the structure all along its life time, it allows to compute the evolution of the reliability in time. The time dependency appearing in these so-called time-variant reliability problems may be of two kinds: • material properties may be decaying in time. Examples of degradation mechanisms are fatigue crack growth in fracture mechanics, corrosion in steel structures or in reinforced concrete rebars, concrete shrinkage and creep phenomena, etc. • loading may be randomly varying in time: in this case, stochastic processes are introduced in the analysis. This allows to account for environmental loading such as wind velocity, temperature, etc. Time-variant reliability aims at computing the evolution in time of probabilities of failure by taking into account the complete evolution of the material

properties, environmental conditions and loading in time. When only degradation mechanisms are considered, the randomness is often represented by random variables multiplied by deterministic functions of time (the “shape” of the degradation kinetics). As these deterioration processes are usually monotonous in time (the crack length increases in time, the size of uncorroded zones decreases in time, etc.), it is possible to focus onto the end of the lifetime of the structure. In this case the problem may be solved using time-invariant reliability tools. The second kind of time dependency has been addressed by many authors in the past two decades, especially Breitung (1988), Schall et al. (1991), Rackwitz (1998) using the so-called outcrossing approach. The probability of failure is related to the mean number of outcrossings of the random process through the limit state surface. The approach, referred to as asymptotic method (AsM) in the sequel, couples analytical results for outcrossings of random processes with asymptotic integration. Hagen and Tvedt (1991) suggested another approach based on a parallel system reliability formulation for computing the outcrossing rate (as defined in the next section). Li and Der Kiureghian (1995), Der Kiureghian (1996) and Der Kiureghian (2000) applied it to the study of a variety of reliability problems such as nonlinear dynamics and inelastic space-

variant finite element analysis. Vijalapura et al. (2000) applied it to the study of the Bouc-Wen single degree-of-freedom hysteretic oscillator. The PHI2 method presented in this paper (Andrieu et al. 2002a, 2002b, Andrieu-Renaud 2002, Sudret et al. 2002) is inspired by this approach. It allows to address in a similar way all kinds of timedependent problems, whatever the time dependency mentioned above. After recalling the basics of time-variant reliability, the PHI2 method is presented. It is then compared to the classical asymptotic approach as implemented in the software COMREL-TV V7.1 (RCP,1998) on two academic examples (stationary and non-stationary) for which analytical solutions can be derived. The quality of the upper bound obtained by the PHI2 method is assessed by Monte Carlo simulation. The numerical features and the computational cost of each method are compared. The PHI2 method is finally applied to study the failure of a concrete beam submitted to creep and stochastic pinpoint load. 2 TIME-VARIANT RELIABILITY CONCEPTS

2.2 The outcrossing approach The outcrossing approach is based on the computation of the outcrossing rate of the process under consideration through the limit state surface. This outcrossing rate is defined by: ν + (τ ) =

lim

∆τ →0 ,∆τ >0

Prob({G (t , X(t ,ω )) > 0}∩ {G (t + ∆τ , X(t + ∆τ ,ω )) > 0}) ∆τ

(4)

This definition requires that the process be regular, that is to say the probability of having more than one crossing in a small interval of time must be negligible. The mean number of crossings in the time interval time [t1,t2] is given then by: (5) t2 E  N ( t 1 , t2 )  = ∫ ν + (τ ) dτ t1 An upper bound to the probability of failure is given by Bolotin (1981): (6) Pf ,c (t1 , t 2 ) ≤ Pf ,i (t1 ) + E [N (t 1, t 2 )] A lower bound is also available (Shinozuka 1965): (7) Pf ,c ( t1 , t2 ) ≥ max Pf ,i (τ ) t1 ≤τ ≤t2

2.1 Problem statement and notations Let us denote by X(t,ω) the set of random variables Xj(ω), j=1,…, p and one-dimensional random processes Xk(t,ω), k=p+1,…, p+q describing the randomness in the mechanical problem under consideration. In this notation, ω stands for the outcome in the space of outcomes Ω. The time-dependent limit state function associated with the reliability analysis is denoted by G(t, X(t,ω)). Positive values of G(t, X(t,ω)) correspond to the safe domain and negative or zero values to the failure domain. Zero values of G(t, X(t,ω)) define the limit state surface. The failure of the structure within the time interval [0,t] is represented by the event: (1) E = {∃τ ∈ [0, t ] G (τ , X(τ ,ω )) ≤ 0}

Thus the probability of failure Pf,c of the structure within the interval [0,t] is defined by: (2) Pf ,c (0, t ) = P(E ) = Prob(∃τ ∈ [0, t ] G (τ , X(τ , ω )) ≤ 0 )

Let us define also the so-called point-in-time (or instantaneous) probability of failure at time t by: (3) Pf ,i ( t ) = Prob G ( t , X ( t , ω ) ) ≤ 0

(

)

This quantity is computed by considering time as a fixed parameter in Eq.(3) and by replacing the random processes by the corresponding random variables.

3 COMPUTING THE OUTCROSSING RATE 3.1 Asymptotic integration Reliability analysis usually involves both S-processes and R-random variables according to the classification by Schall et al. (1991). Thus the results mentioned above are conditioned to values of {Rj, j = 1,…,p}. The upper bound to the probability of failure (6) is computed by integrating the latter both over all realizations of the R’s and over time. Breitung (1988), Bryla et al. (1991), Hagen (1992) used this theory and asymptotic integration to solve the problem for scalar Gaussian processes using Laplace integration and FORM-like procedures. COMREL (RCP 1998) implements this kind of approach and generalizes it to vector processes. 3.2 The PHI2 method The outcrossing rate is computed by the twocomponent parallel system analysis (Eq.(4)), i.e. two successive time-invariant analysis (using simulation or FORM) and the evaluation of the binormal cumulative distribution function as presented first by Hagen et Tvedt (1991) and then used by Li and Der Kiureghian (1995). The difference with the previous applications of this approach is that we consider that the process involved in the study is the limit-state function itself. We do not make a first order Taylor expansion of the limit-state at t+∆t (See Eq.(4)).

Moreover the random process under consideration should not be discretized as in Der Kiureghian (2000). Let us denote by: • β(t) the (time-invariant) reliability index associated with the event {G(t, X(t,ω))>0}, and α (t ) the unit normal to the limit state surface at the design point. • β(t+∆t) the (time-invariant) reliability index associated with the event {G(t+∆t, X(t+∆t,ω))≤0} and α (t + ∆t ) the unit normal to the limit state surface at the design point. • ρG (t , t + ∆t ) = −α (t ) ⋅α (t + ∆t ) the correlation between the two events { G(t, X(t,ω))>0} and { G(t+∆t, X(t+∆t,ω))≤0} The first order evaluation of the probability of failure of the parallel system writes: (8) Φ (β (t ),− β (t + ∆t ); ρ G (t , t + ∆t )) ν (τ ) = 2 ∆t

The choice of ∆τ is crucial. It is discussed in (Andrieu et al. 2002b). In case of stationary problems, the upper bound to Pf is obtained by multiplying ν by the time interval under consideration. When non-stationarity is present, the outcrossing rate is computed at several time instants within the time interval and the integration in Eq.(5) is carried out using Simpson’s rule. 3.3 Monte-Carlo simulation Monte-Carlo simulation is a general method to solve reliability problems that gives accurate results if the number of samples is large enough. In case when a random process appears in the limit state function, a strategy to generate trajectories of the process has to be selected. The usual approach is based on the random process discretization, i.e. its representation by a finite set of random variables. Several schemes are presented in Schüeller et al. (1997) and Sudret & Der Kiureghian (2000). The method retained in the present paper is the so-called Expansion Optimal Linear Estimation method (EOLE) (Li & Der Kiureghian 1993). It is computed by a MATLAB routine presented in Sudret & Der Kiureghian (2000), which can be http://www.ce.berkeley.edu/ downloaded at ~haukaas under the item FERUM. An error estimator εrr(t) which allows to evaluate the accuracy of the discretization is also given in Sudret & Der Kiureghian (2000). Suppose a trajectory of the limit state function G(t, x(t,ω)) is simulated in the time interval [0,t]. Its values are stored in an array {G(ti), ti=(i-1).∆t, i=1,…,N}, where N is the number of points used in the discretization and ∆t=t/(N-1) is the sampling step. Let ki, i=1,…,N-1 be a set of counters that are incremented by one each time the current trajectory presents an outcrossing of the limit state surface

within the interval [ti,ti+1]. The estimated probability of failure is given by: Pf ,c (0, ti ) =

1 N MC

i −1

∑k

j

, 2≤i≤ N

(9)

j =1

The variance of this estimator is: Var (Pf ,c (0, t i )) =

1 N MC

 1  − 1  N MC 

 1 k j −  ∑ j =1  N MC i −1

 k j  ∑ j =1  i −1

2

  

(10)

It will be used in the sequel to evaluate the accuracy of the simulations. 4 BENCHMARK OF THE APPROACHES In this section, two examples are dealt with to benchmark various approaches: • the Asymptotic Method (AsM), as implemented in COMREL-TV V.7.1 (RCP 1998), • the PHI2 method. COMREL-TI is used to carry out the required time-invariant FORM analyses. Note that the default parameters for FORM analysis are used, except the tolerance on the convergence, which is set equal to 10-3. • Monte Carlo simulation: MATLAB and MathCad routines are used to discretize the random process and perform the simulation. 4.1 Example #1: stationary case R(ω)-S(t,ω) The first example deals with the outcrossing of a Gaussian random threshold R(ω) (mean µR, standard deviation σR) by a Gaussian stationary random process S(t,ω) (mean µS, standard deviation σS). The parameters defining R and S are given in Table 1, column #2. Table 1: Definition of variables R and S for the various case studies Variable Example #1 Example #2 R(ω) µR=150 µR=150 σR=15 σR=15 S(t,ω) µS=50 µS=5 σS=17.5 σS=17.5 a -2.5

The autocorrelation coefficient function of the process is denoted by ρS(t1,t2): (11)  (t 2 − t1 )2   ρ S (t1 , t 2 ) = exp − l 2   where the correlation length l of the process is set equal either to 0.5 or 1. The corresponding cycle rate is ω 0 = 2 / l . The time interval [0,T] under consideration is [0,30]. The limit state function under consideration is: (12) G(t , X(t , ω )) = R(ω ) − S (t , ω )

As the case is stationary, only the outcrossing rate ν+ will be studied. The closed-form solution for ν+ was derived in Sudret et al. (2002) based on Rice’s formula:

length is 10 m, its cross-section is rectangular (b = 0.2 m, h = 0.2 m). The limit stress is fc,28 = 30 MPa. The concrete instantaneous Young’s modulus is E28 = 11000 (fc,28)1/3 = 34180 MPa. This beam is subjected to dead load (denoting by ρ = 2500 kg/m3 the concrete mass density, this load is equal to p = ρbh=1000 N/m) as well as a pinpoint load S applied onto its middle point. The maximum deflection in this point due to the dead load is: (16) 5 p l4 0 δw = 384 E28 I where I=bh3/12 is the beam inertia. Let us denote by J(t,tc) the uniaxial compliance of concrete, which depends on time t as well as time tc corresponding to the loading of the beam. According to linear visco-elasticity, the uniaxial strain can be computed as:

ε ( t ) = J ( t , tc ) .σ

(17)

The French codified model for creep is cast as (BPEL, 1991): J ( t , tc ) =

1 (1 + ϕ ( t , tc ) ) E28

(18)

where ϕ ( t , tc ) is the creep coefficient. Using this notation, the evolution in time of the maximal deflection is: 5 p l4 δ w (t ) = J (t , t c ) 384 I =δ

0 w

(19)

(1 + ϕ ( t , t ) ) c

The creep model used in the sequel is defined by:

ϕ ( t , tc ) = ϕ ∞

t − tc

(20)

t − tc + K

The values of the different parameters are: •

tc=28 days,

•

K=15.8 is computed from the characteristics of the cross-section of the beam.

•

ϕ∞ = 2

The live load is a pinpoint load S applied at midspan during one year starting from time t1. The resulting deflection at this point is: Sl 3 δ S (t ) =  H ( t − t1 ) − H ( t − t2 )  48 E28 I 

(21)

where H is t2 = t1 + 1 year .

the

Heaviside

function

and

5.2 Reliability problem statement Let us denote by δ lim the maximal acceptable deflection. The limit state function associated with the exceed of δ lim reads: G ( t , X ( t , ω ) ) = δ lim − δ w ( t ) − δ S ( t )

(22)

In this study, we suppose that the creep phenomenon (ϕ∞ , K) as well as the loading (p, S, t1) are random. The input parameters are gathered in Table 5. Table 5: Random variables and parameters for the concrete beam under creep and pinpoint load Parameters Type of law Mean c.o.v p Lognormal 1000 N/m 10 % ϕ∞ Lognormal 2 20 % K Lognormal 15.8 30 % S Lognormal 5000 N 20 % t1 Uniform [0,49 years]

The study is carried out over [0,50 years]. 5.3 Details of the calculation The following steps are necessary: 1. compute the reliability index β(t) and β(t+∆t) at different instants t and t+∆t 2. evaluate the correlation ρG(t,t+∆t) between the limit state at t and that at t + ∆t . 3. compute the outcrossing rate and then the failure probability by integration. 5.3.1 Computing the reliability index β(t) and β(t+∆t) For practical computation, the limit state function is rewritten as:  δ − δ w ( t ) for t < t1 or t > t2 G ( t , X ( t , ω ) ) =  lim δ lim − δ w ( t ) − δ S ( t ) for t1 ≤ t ≤ t2

(23)

The FORM method could not be applied for computing the point-in-time reliability index β(t). Indeed, during the search of the design point, one or the other expression of (23) is used depending on the current realization of t1, meaning a switch between two different states. Moreover, the gradient ∂G/∂t1 is required and is always equal to zero (precisely it would be a Dirac function which cannot be handled properly). A Monte Carlo simulation strategy has eventually been adopted to get β(t) and β(t+∆t). Note that point-in-time problems are considered here, meaning that the process S is not discretized (in contrast to

length is 10 m, its cross-section is rectangular (b = 0.2 m, h = 0.2 m). The limit stress is fc,28 = 30 MPa. The concrete instantaneous Young’s modulus is E28 = 11000 (fc,28)1/3 = 34180 MPa. This beam is subjected to dead load (denoting by ρ = 2500 kg/m3 the concrete mass density, this load is equal to p = ρbh=1000 N/m) as well as a pinpoint load S applied onto its middle point. The maximum deflection in this point due to the dead load is: (16) 5 p l4 0 δw = 384 E28 I where I=bh3/12 is the beam inertia. Let us denote by J(t,tc) the uniaxial compliance of concrete, which depends on time t as well as time tc corresponding to the loading of the beam. According to linear visco-elasticity, the uniaxial strain can be computed as:

ε ( t ) = J ( t , tc ) .σ

(17)

The French codified model for creep is cast as (BPEL, 1991): J ( t , tc ) =

1 (1 + ϕ ( t , tc ) ) E28

(18)

where ϕ ( t , tc ) is the creep coefficient. Using this notation, the evolution in time of the maximal deflection is: 5 p l4 δ w (t ) = J (t , t c ) 384 I =δ

0 w

(19)

(1 + ϕ ( t , t ) ) c

The creep model used in the sequel is defined by:

ϕ ( t , tc ) = ϕ ∞

t − tc

(20)

t − tc + K

The values of the different parameters are: •

tc=28 days,

•

K=15.8 is computed from the characteristics of the cross-section of the beam.

•

ϕ∞ = 2

The live load is a pinpoint load S applied at midspan during one year starting from time t1. The resulting deflection at this point is: Sl 3 δ S (t ) =  H ( t − t1 ) − H ( t − t2 )  48 E28 I 

(21)

where H is t2 = t1 + 1 year .

the

Heaviside

function

and

5.2 Reliability problem statement Let us denote by δ lim the maximal acceptable deflection. The limit state function associated with the exceed of δ lim reads: G ( t , X ( t , ω ) ) = δ lim − δ w ( t ) − δ S ( t )

(22)

In this study, we suppose that the creep phenomenon (ϕ∞ , K) as well as the loading (p, S, t1) are random. The input parameters are gathered in Table 5. Table 5: Random variables and parameters for the concrete beam under creep and pinpoint load Parameters Type of law Mean c.o.v p Lognormal 1000 N/m 10 % ϕ∞ Lognormal 2 20 % K Lognormal 15.8 30 % S Lognormal 5000 N 20 % t1 Uniform [0,49 years]

The study is carried out over [0,50 years]. 5.3 Details of the calculation The following steps are necessary: 1. compute the reliability index β(t) and β(t+∆t) at different instants t and t+∆t 2. evaluate the correlation ρG(t,t+∆t) between the limit state at t and that at t + ∆t . 3. compute the outcrossing rate and then the failure probability by integration. 5.3.1 Computing the reliability index β(t) and β(t+∆t) For practical computation, the limit state function is rewritten as:  δ − δ w ( t ) for t < t1 or t > t2 G ( t , X ( t , ω ) ) =  lim δ lim − δ w ( t ) − δ S ( t ) for t1 ≤ t ≤ t2

(23)

The FORM method could not be applied for computing the point-in-time reliability index β(t). Indeed, during the search of the design point, one or the other expression of (23) is used depending on the current realization of t1, meaning a switch between two different states. Moreover, the gradient ∂G/∂t1 is required and is always equal to zero (precisely it would be a Dirac function which cannot be handled properly). A Monte Carlo simulation strategy has eventually been adopted to get β(t) and β(t+∆t). Note that point-in-time problems are considered here, meaning that the process S is not discretized (in contrast to

the Monte Carlo simulation of the time-variant problems presented in section 4) but replaced by appropriate random variables. Precisely, let us denote by: •

St the random variable associated to the process at time t • St+∆t the random variable associated to the process at time t+∆t, not correlated with St • S~t + ∆ t the random variable associated to the process at time t+∆t, correlated with St by the autocorrelation coefficient function ρ S . The following relationship holds: (24) St +∆t = ρ S ( t , t + ∆t ) St + 1 − ρ S 2 ( t , t + ∆t ) St +∆t

The latter equation is used to simulate properly correlated realizations of ( St , St +∆ t ) from the generated independent samples ( St , St +∆ t ) .

A total number of 1,000,000 samples have been used at each time step t. This high number was necessary to obtain a coefficient of variation (c.o.v) less than 5%. The obtained point-in-time reliability index is plotted as a function of time in Figure 3, the coefficient of variation of the simulation being plotted in Figure 4. 6 5 4

β (t)

3 2 1 0 0

10

20

30

40

50

Time (years)

Figure 3: Concrete beam submitted to creep and pinpoint load: point-in-time reliability index vs. time

5.3.2 Computation of the correlation ρG(t, t+∆t) As FORM could not be used, the correlation ρG(t, t+∆t) used in Eq.(8) could not be obtained from the unit normal vectors α (t ) . However it can be obtained by carrying out a linear regression analysis (Saporta, 1990) of samples

{G (t + ∆t, x

(i )

(t + ∆t )

)}

vs. samples

{G (t, x (t ) )} . (i )

Figure 5 illustrates the method at

t = 10 years.