Evaluation of the Accident Rate of a Plant Equipped with an Aging Single Protective Channel by the Method of Supplementary Variables

ABSTRACT

This work calculates the reliability of protective systems of industrial facilities, such as nuclear, to analyze the case of equipment subject to aging, important in the extension of the qualified life of the facilities. By means of the method of supplementary variables, a system of partial and ordinary integral-differential equations was developed for the probabilities of a protective system of an aging channel. The system of equations was solved by finite differences. The method was validated by comparison with channel results with exponential failure times. The method of supplementary variables exhibits reasonable results for values of reliability attributes typical of industrial facilities.

KEYWORDS:
industrial facilities; protective systems; non-markovian processes; Markov chains; finite differences; useful life extension

1 INTRODUCTION

The reliability analysis of protective system components is of paramount importance and widely studied. This relevance is related to the fact that the protection function in an industrial installation is fundamental for its integrity, because it monitors process parameters, such as pressures and temperatures and, if necessary, commands the shutdown of the facility. Classically, this problem has been addressed by considering its behavior when the plant it should protect is subject to high demand rates⁶ and, later, the influence of the channel repair rate was addressed⁹ and also by considering a two redundant channel protective system¹⁰. A further extension up to 5 redundant protective channels has been considered⁵, where different actuation (voting) logics have been taken into account.

However, it is necessary to address the components that are no longer in their useful life, that is, they are aging and this is the purpose of this paper. When a protective channel is not in its useful life anymore, two options are normally considered: 1) replace the aged channel; 2) analyze the possible extension of the aging channel life. Next, a cost-benefit analysis is performed and the proper decision is made.

For protection systems, the reliability parameter of interest is its average unavailability, U, which depends on the failure, λ, and repair rates, µ, on the number of channels constituting it, in addition to the demand rate, ν, and also on the test and maintenance policies and their actuation logic. However, the characteristic considered is the plant accident rate, η⁶, which is defined as

$\eta =\nu U(\lambda ,\mu ,\nu ).$

For an initial evaluation, we consider a protective system made by a single aging channel. Figure 1 shows the diagram that represents the transition logic between the possible system states. P _i (x) represents the probability that the channel is in the i −th state, x is the channel age; λ(x) stands for the age-dependent failure rate; µ is the channel repair rate; ν is the channel demand rate; and γ is the probability of a human error during the channel maintenance. This probability models the situation in which the channel repair can induce its failure due to human error. The indices in the triplet ¡i, j, k¿ represent the number of channels in operation, the number of failed channels with undisclosed failures and the number of failed channels detected, respectively. The distinction between covered and uncovered faults is important in this context, since channel failure detection only occurs in two situations: when the system is called for (a real situation), or when the system is tested.

Figure 1
State transition diagram for an aging channel.

Defining τ _p as the channel proof-test interval, the plant accident rate will be given by

$\eta =\frac{\nu }{{\tau }_p}{\int }_0^{{\tau }_p}{P}_2\left(t\right)dt$ (1.1)

where, as defined earlier, P ₂(t) represents the probability of the system meeting the faulty channel (i.e., in state 2) and its failure has not been detected. Eq. (1.1) means that the only possible accident initiating events are those caused by demands of the system while it is faulty and its failure has not been detected, since it is assumed that during repair the plant is off, which means we are considering offline repair only.

2 THE METHOD OF SUPPLEMENTARY VARIABLES

Even if one of the transition rates is a function of time (calendar time t, in general), the model can still be considered Markovian by inserting a supplementary variable (x), which takes into account the channel age. In this case, it is necessary to include it in the analysis, since the absence of the memory property is no longer valid because the channel failure times (t) do not follow exponential distributions anymore. The system of integral-differential equations^{2), (12} will be

$\frac{\partial p_1(x,t)}{\partial x}+\frac{\partial p_1(x,t)}{\partial t}=-\lambda \left(x\right)p_1(x,t),$ (2.1)

$\frac{d{P}_2\left(t\right)}{dt}={\int }_0^{\infty }\lambda \left(x\right)p_1(x,t)dx-\nu P_2\left(t\right)+\gamma \mu P_3\left(t\right),$ (2.2)

$\frac{d{P}_3\left(t\right)}{dt}=\nu P_2\left(t\right)-\mu P_3\left(t\right).$ (2.3)

The states presented in Figure (1) have the following meanings: state 1 (system in operation) is modeled by means of a probability density $p_1(x,t)$, interpreted as: $p_1(x,t)\Delta x\Delta t$ is the probability that the channel is in state 1 between the instants of time t and $t+\Delta t$ with age between x and $x+\Delta x$. The probability of the system being in state 1 at time t will be given by

$P_1\left(t\right)={\int }_0^{\infty }{p}_1(x,t)dx.$ (2.4)

State 2 represents the system when the channel is failed and its failure has not been revealed. State 3 represents the failed system, with the failure revealed.

The initial condition for Eq. (2.1) is given by

$p_1(x,0)=f\left(x\right),$ (2.5)

where f(x) is the probability density function of the channel failure times. In this work, we admit that the failure times follow a three-parameter Weibull distribution and, therefore, we will have for λ(x) and f(x), respectively⁷.

$\lambda \left(x\right)=\left\{\begin{array}{cc}{\lambda }_0,& \begin{array}{cc}\text{if}& x<x_0\end{array}\\ {\lambda }_0+\left(\frac{m}{\theta }\right){\left(\frac{x-x_0}{\theta }\right)}^{m-1},& \text{if}x\ge x_0\end{array}\right.$ (2.6)

and

$f\left(x\right)=\left\{\begin{array}{cc}{\lambda }_0{e}^{-{\lambda }_{0}x},& \text{if}x<x_0\\ \begin{array}{c}\left[{\lambda }_0+\left(\frac{m}{\theta }\right){\left(\frac{x-x_0}{\theta }\right)}^{m-1}\right]e^{-{\lambda }_{0}x-{\left(\frac{x-x_0}{\theta }\right)}^m},\end{array}& \text{if}x\ge x_0\end{array},\right.$ (2.7)

where m is the shape parameter and $m>1$, as the system ages. θ represents the scale parameter of the Weibull distribution; λ ₀ is the failure rate of the system before it ages, and therefore is constant and x ₀, the age of the channel at the beginning of its aging period. The Weibull distribution is commonly used for representing failure times under aging¹². Figure (2) shows the behavior of the channel failure rate, Eq. (2.6).

Figure 2
Behavior of the channel failure rate wih time.

It is clear from Figure (2) that the channel failure rate is equal to $1/yr.$ (i.e., constant) for the first year of the analysis $(x_0=1yr.)$ and then, presents an increasing concave upward behavior for the second year of the analysis. Notice that we are not considering in our analysis the wear-in period of the bathtub curve¹², just the useful life period (constant failure rate model) and the wear-out period (monotonically increasing failure rate). We assume that the channel useful life period starts at $t=0$.

Figure (3) shows the behavior of the channel failure density, given by Eq. (2.7). Since it represents the useful life region and, from a given point (represented by the channel age) the behavior in the aging period (which, in the case of the figure, was considered equal to 1 year), there is a small change in the derivative of this function. Besides Eq. (2.5), we also need the following conditions

$p_1(x_0,t)=(1-\gamma )\mu P_3\left(t\right),P_1\left(0\right)=1,P_i\left(0\right)=0,\mathrm{i}=2,3.$ (2.8)

Figure 3
Behavior of the channel failure density.

Eq. (2.8) displays three important information. The boundary condition shows the state of the channel at age x ₀. Probabilities $P_i\left(x_0\right)\ne Pi\left(0\right)(i=1,2,3)$ are known, because before aging the failure times follow an exponential distribution. The second and third information are the probabilities of the protection system to be found in each state, 1, 2 or 3, at the beginning of its operation, that is, at $t=0$. This means that at this time the system is working. Of course, $\sum _{i=1}^3{P}_i\left(t\right)=1$. From the definition of Eq. (1.1), we note that we are interested only in the probability P ₂(t), since, as already mentioned, we consider that only offline repairs are performed.

3 SOLUTION OF THE SYSTEM OF EQUATIONS

The system of equations (2.1)-(2.5) and (2.8) is coupled. Thus, the solution method chosen was that of finite differences¹. For the numerical solution of the integrals in the equations of the problem, the Simpson Composite or 1/3 method of Simpson was used¹.

The first step is to define the mesh over which the method will work

$x_1=x_0,x_{j+1}=x_0+j\Delta x,x_{j+1}-x_j=\Delta x;\forall j=1,\cdots ,J$

and

$t_1=0,t_{l+1}=l\Delta t,t_{l+1}-t_l=\Delta t;\forall l=1,\cdots ,L,$

where index j refers to the supplementary variable (channel age), J is the number of steps considered for the channel age interval and index l represents calendar time, where L is the number of steps considered for the calendar time interval. Extending this notation for probabilities, we will have

$p_1(x_j,t_l)=p_{1,j}^l,\forall j=1,\cdots ,J+1;l=1,\cdots ,L+1$

and

$P_i\left(t_l\right)=P_i^l,l=1,\cdots ,L+1,i=2,3.$

Applying the boundary conditions, Eqs. (2.5) and (2.8), we shall have, respectively

$$\begin{gathered} p_{1,j}^1=f_j,\forall j=1,\cdots ,J+1, \\ {p}_{1,1}^l=(1-\gamma )\mu P_3^l,\forall l=2,\cdots ,L+1, \\ {P}_1^1=1, \end{gathered}$$

and

$P_i^1=0,i=2,3.$

For the discretization of Eq. (2.1), we will make the following approximation³

$p_{1,j}^l\simeq \frac{1}{2}\left[p_{1,j+1}^l+p_{1,j-1}^l\right].$ (3.1)

We will have, for the partial derivatives

$$\begin{gathered} \frac{\partial p_1(x,t)}{\partial x}\simeq \frac{p_{1,j+1}^l-p_{1,j-1}^l}{2\Delta x}, \\ \frac{\partial p_1(x,t)}{\partial t}\simeq \frac{p_{1,j}^{l+1}-p_{1,j}^l}{\Delta t}. \end{gathered}$$ (3.2)

Thus, putting Eq. (3.1) into Eq. (3.2), we will have

$\frac{\partial p_1(x,t)}{\partial t}\simeq \frac{p_{1,j}^{l+1}-\frac{1}{2}\left[p_{1,j+1}^l+p_{1,j-1}^l\right]}{\Delta t}.$ (3.3)

This discretization is known as the Lax-Friedrichs method and has conditional stability³.

Thus, Eq. (2.1) will become

$p_{1,j}^{l+1}=-{\lambda }_j\Delta t{p}_{1,j}^l+\frac{1}{2}\left(1+\frac{\Delta t}{\Delta x}\right)p_{1,j-1}^l+\frac{1}{2}\left(1-\frac{\Delta t}{\Delta x}\right)p_{1,j+1}^l.$ (3.4)

This equation is valid for $j=2,\cdots ,J;l=1,\cdots ,L$. On the other hand, for $j=J+1$, we use Eq. (3.1) in the following way

$p_{1,J+1}^{l+1}=2p_{1,J}^{l+1}-p_{1,J-1}^{l+1}.$ (3.5)

We will then have, for the derivatives

$\frac{d{P}_i\left(t\right)}{dt}\simeq \frac{P_i^{l+1}-P_i^l}{\Delta t};\forall l=1,\cdots ,L,i=2,3.$ (3.6)

In Eq. (2.2) we have to approximate the integral numerically. As previously mentioned, the method chosen was that of coumpound Simpson¹, so we will have

$$\begin{gathered} {\int }_0^{\infty }\lambda \left(x\right)p_1\left(x,t\right)dx\simeq \\ \frac{\Delta x}{3}\left({\lambda }_1{p}_{1,1}^l+4\sum _{j=1}^{\frac{J}{2}}{\lambda }_{2j}{p}_{1,2j}^l+2\sum _{j=1}^{\frac{J}{2}-1}{\lambda }_{2j+1}{p}_{1,2j+1}^l+{\lambda }_{J+1}{p}_{1,J+1}^l\right). \end{gathered}$$ (3.7)

for even J. Putting Eq.(3.7) and (3.6) into Eq.(2.2), we will have

$P_2^{l+1}=P_2^l+\Delta t\left(-\nu P_2^l+\gamma \mu P_3^l+I\right)$ (3.8)

where

$I=\frac{{\displaystyle \Delta x}}{{\displaystyle 3}}\left({\lambda }_1{p}_{1,1}^l+4\sum _{j=1}^{\frac{{\displaystyle J}}{{\displaystyle 2}}}{\lambda }_{2j}{p}_{1,2j}^l+2\sum _{j=1}^{\frac{{\displaystyle J}}{{\displaystyle 2}}-1}{\lambda }_{2j+1}{p}_{1,2j+1}^l+{\lambda }_{J+1}{p}_{1,J+1}^l\right).$

Finally, for probability P ₃(t), putting Eq. (3.6) into Eq. (2.3), we will have

$P_3^{l+1}=P_3^l+\Delta t\left(-\mu P_3^l+\nu P_2^l\right).$ (3.9)

In order to determine P ₁(t), we apply the Simpson compound method to Eq (2.4)

$P_1^l=\frac{{\displaystyle \Delta x}}{{\displaystyle 3}}\left(p_{1,1}^l+4\sum _{j=1}^{\frac{{\displaystyle J}}{{\displaystyle 2}}}{p}_{1,2j}^l+2\sum _{j=1}^{\frac{{\displaystyle J}}{{\displaystyle 2}}-1}{p}_{1,2j+1}^l+p_{1,J+1}^l\right).$

By putting $l=1$ in Eq. (2.1) we have

$$\begin{gathered} p_{1,1}^1=f_1,p_{1,1}^2=(1-\gamma )\mu P_3^2, \\ {p}_{1,j}^2=-{\lambda }_j\Delta t{p}_{1,j}^1+\frac{{\displaystyle 1}}{{\displaystyle 2}}\left(1+\frac{{\displaystyle \Delta t}}{{\displaystyle \Delta x}}\right)p_{1,j-1}^1+\frac{{\displaystyle 1}}{{\displaystyle 2}}\left(1-\frac{{\displaystyle \Delta t}}{{\displaystyle \Delta x}}\right)p_{1,j+1}^1,\forall j=2,\cdots ,J \end{gathered}$$

and

$p_{1,J+1}^2=2p_{1,J}^2-p_{1,J}^2,j=J+1.$

Depending on the values of ∆t and ∆x adopted, the solution may become unstable. For this system of equations, the method is stable if

$\frac{{\displaystyle \Delta t}}{{\displaystyle \Delta x}}\le 1,$ (3.10)

², obtained through the analysis of Eq. (3.4). Besides this condition, we must also have: $\Delta t\le 1/\nu$ and $\Delta t\le 1/\mu$. These conditions come from the analysis of Eqs. Eq.(3.8) and Eq.(3.9). Therefore, in the first place, the value of ∆t is chosen from:

$\Delta t=\text{min}\left(\frac{{\displaystyle 1}}{{\displaystyle \nu }},\frac{{\displaystyle 1}}{{\displaystyle \mu }}\right)$ (3.11)

and next, we use Eq. (3.10) in order to choice the channel age step, ∆x. With the solutions to the probabilities, we can then calculate the plant accident rate, again using the method of compound Simpson:

$\eta =\frac{{\displaystyle \nu }}{{\displaystyle {\tau }_p}}\frac{{\displaystyle \Delta t}}{{\displaystyle 3}}\left(P_2^1+4\sum _{l=1}^{\frac{{\displaystyle L}}{{\displaystyle 2}}}{P}_2^{2l}+2\sum _{l=1}^{\frac{{\displaystyle L}}{{\displaystyle 2}}-1}{P}_2^{2l+1}+P_2^{L+1}\right).$ (3.12)

4 RESULTS

Table 1 displays the input data used for the accident frequency calculations. Those data that are without a bibliographic reference were chosen based on our experience. As we are dealing with monitoring and shutdown (thus, protective) systems, the search for failure data is much more complicated since they can be applied to any kind of industrial plant. In this sense, we chose some of the input data displayed in Table 1 based on our field experience without discriminating from what kind of plant they have come from. Besides, in many cases these data are proprietary information.

Table 1
Values used for the calculations of the reliability figures of merit.

Parameter	Value	Parameter	Value
ν	10 yr.−1(1	θ	1 yr.
µ	52 yr.−1(1	∆t	5 E-05 yr.
γ	0.1	∆x	5 E-05 yr.
τ p	1 yr.1	x 0	1 yr.
λ 0	1.0 yr.−1(1	t L+1	2 yr.
m	2.5	x J+1	50 yr.

Figure (4) shows the behavior of P ₂(t), the probability of the channel presenting an unrevealed failure. This probability is important for calculating the plant accident rate of the facility, according to Eq. (3.12). It is seen from this figure that the probability reaches a steady value around one year (using the data presented in Table 1). An important aspect in this context is that the interval between channel tests is not greater than one year in practice. Therefore, to approximate the probability P ₂ by its asymptotic value (which simplifies considerably the calculations), would generate a pessimistic result (that is, a frequency of accidents greater than the real one). In this sense, it is not generally advisable to work with asymptotic probabilities in this context.

Figure 4
Behavior of the channel probability of undergoing an unrevealed failure over time.

We have tested this method earlier for exponential times¹¹.

The results of the analysis are presented in Table 2. Given the behavior of the failure rate [Eq. (2.6)] the results for the frequency of accidents are larger than the corresponding ones for the useful life of the plant. This is expected behavior. Calculating the frequency by the supplementary variable method shows that it is possible to stipulate realistic inter-aging periods and thus establish a system useful life, which has important cost implications. The parameter θ ( channel characteristic life¹²) is directly proportional to the channel mean time to failure. It is clear from Table 2 that the lesser the characteristic channel life, the greater the plant accident rate, as expected.

Table 2
Accident rates for the plant with an aging channel

θ (yr.)	ν (yr. −1 )	η(yr. −1 )
10	0.5	0.07
	10.0	0.21
	100.0	0.22
1.0	0.5	0.23
	10.0	0.93
	100.0	1.07

Fig. 5 shows the behavior of the accident rate for two scenarios: in the first, the characteristic life of the channel is equal to 1 year and, in the second one, to 10 years.

Figure 5
Plant accident rate with an aging single channel.

In the first case, it is noted that the accident rate reaches a stationary value for a demand rate of about 70/year, while in the second one, for a demand rate of about 5/yr.

The investigated solution method returns reasonable values for demand rates of all ranges. Considering lower demand rates means making use of very high demand times since these follow an exponential distribution. According to the number of points with which the mesh is defined, J and L, the computation time greatly increases, with L being the most affected. For <mml:math><mml:mi>L</mml:mi><mml:mo> </mml:mo><mml:mo>></mml:mo><mml:mo> </mml:mo><mml:mn>106</mml:mn></mml:math>, the calculation time is of the order of 10 h. This entails two needs: that of considering the mesh refinement, until an optimal one is found, and that of keeping the interval of time t short. The optimization of the mesh size was established so that the result of the integral of the probability density function had an error smaller than 10⁻¹⁰ . This error value was chosen because of the precision of the method. It was concluded that for errors of order of magnitude greater than 10⁻⁶ , the method becomes rapidly unstable. The choice of errors of lesser order of magnitude does not influence the calculation time, since the integrals that demand greater computational effort are all in relation to x, that is, they depend on J.

5 CONCLUSIONS

The purpose of this paper was to evaluate the accident rate of a plant equipped with a singlechannel protective system when this protective channel is under aging. The developed model can serve as a useful tool to assist in the decision making of repair policies and life extension of an industrial facility in a cost-efficient way. The results presented show a behavior compatible with that of the channel in its aging period. The possibility of using stationary values of the channel undetected probability of failure (this means the trip channel is in state 2 of Fig. 1) has shown that, although the calculations involved are simpler, they produce conservative (i.e., pessimistic) accident rates. These conservative results have an important impact over the repair policies to be set for the plant, in the sense that they can burden the number of man-hours expended. It is intended to analyze the trip system with redundancy, because it is a common practice in industrial facilities: a redundancy can reduce the chance of plant shutdown failure when needed. Application to systems with up to 5 channels covers 90% of the protection systems used in practice^{5), (10}. In this case, possible redundancy alternatives can be evaluated in comparison to the same configurations over their lifetime. Comparison of the method applied here with others, for example that of the stages^{2), (12}, is the subject of subsequent research. The application of the Generalized Perturbation Theory^{4), (8} would provide a better understanding of the system behavior and its peculiarities, allowing the development of sensitivity analyses of the parameters involved. This feature is quite important because data variability is to be expected in real cases and the plant accident rate behavior is an important piece of information for establishing maintenance policies.

ACKNOWLEDGMENTS

We would like to thank our colleague Dr. M. C. Martins for his invaluable help with the figures.

REFERENCES

• ¹
  A. Alvim. “Mumerical Methods in Nuclear Engineering (in Portuguese)”. Ed. Certa, Curitiba (2007).
• ²
  D.R. Cox & D.H. Miller. “The Theory of Stochastic Processes”. Wiley, London (1965).
• ³
  A. Fortuna. “Computational Techniques for Fluid Dynamics: Basic Concepts and Applications (in Portuguese)”. Edusp, São Paulo (2000).
• ⁴
  P.F. Frutuoso e Melo, A.C.M. Alvim & F.C. Silva. Sensitivity Analysis on the Accident Rate of a Plant Equipped with a Single Protective Channel by Generalized Perturbation Methods. An. Nucl. Energy, 25 (1998), 1191-1207.
• ⁵
  P.F. Frutuoso e Melo, L.F. Oliveira & R. Youngblood. A Markovian model for the reliability analysis of multi-channeled protective systems considering revealed failures and common-cause failures by the alpha model. Proc. 9th Nat. Meet. On Reactor Physics and Thermal Hydraulics, Braz. Assoc. Nucl. Energy, (1993), 440-446.
• ⁶
  F.P. Lees. A General Relation for the Reliability of a Single-Channel Trip System. Reliability Engineering, 3(1) (1982), 1.
• ⁷
  E.E. Lewis. “Introduction to Reliability Engineering”. New York: Wiley, New York (1996).
• ⁸
  E.F. Lima, D.G. Teixeira, P.F. Frutuoso e Melo, F.C. Silva & A.C.M. Alvim. Sensitivity Analysis of the Accident Rate of a Plant by the Generalized Perturbation Theory. International Journal of Mathematical Models and Methods in Applied Sciences, 10 (2016), 309-316.
• ⁹
  L.F. Oliveira & J.D.S. Neto. Influence of the Demand Rate and Repair Rate on the Reliability of a Single-Channel Protective System. Reliability Engineering, 17 (1987), 267-276.
• ¹⁰
  L.F. Oliveira, R.W. Youngblood & P.F. Frutuoso e Melo. Hazard Rate of a Plant Equipped with a Two-Channel Protective System Subject to a High Demand Rate. Reliability Engineering and System Safety, 28 (1990), 35-58.
• ¹¹
  M.O. Pinho, H.C.N. Fernandez, A.C.M. Alvim & P.F. Frutuoso e Melo. Availability of a Component Subject to an Erlangian Failure Model Under Wearout by Supplementary Variables. J. of the Braz. Soc. Mechanical Sciences, 21(1) (1999), 109-122.
• ¹²
  C. Singh & R. Billinton. “System Reliability Modelling and Evaluation”. Hutchinson, London (1977).

Publication Dates

• Publication in this collection
  08 Nov 2021
• Date of issue
  Oct-Dec 2021

History

• Received
  28 Jan 2021
• Accepted
  16 Mar 2021