Dynamic reliability of digital-based transmitters

Reliability Engineering and System Safety 96 (2011) 793–813

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Dynamic reliability of digital-based transmitters Florent Brissaud a,c,n, Carol Smidts b, Anne Barros c, Christophe Be´renguer c a b c

Institut National de l’Environnement Industriel et des Risques (INERIS), Parc Technologique Alata, BP 2, 60550 Verneuil-en-Halatte, France The Ohio State University (OSU), Nuclear Engineering Program, Department of Mechanical Engineering, Scott Laboratory, 201 W 19th Ave, Columbus OH 43210, USA Universite´ de Technologie de Troyes (UTT), Institut Charles Delaunay (ICD) & UMR CNRS 6279 STMR, 12 rue Marie Curie, BP 2060, 10010 Troyes Cedex, France

a r t i c l e i n f o

abstract

Article history: Received 25 May 2010 Received in revised form 12 December 2010 Accepted 24 December 2010 Available online 1 January 2011

Dynamic reliability explicitly handles the interactions between the stochastic behaviour of system components and the deterministic behaviour of process variables. While dynamic reliability provides a more efficient and realistic way to perform probabilistic risk assessment than ‘‘static’’ approaches, its industrial level applications are still limited. Factors contributing to this situation are the inherent complexity of the theory and the lack of a generic platform. More recently the increased use of digitalbased systems has also introduced additional modelling challenges related to specific interactions between system components. Typical examples are the ‘‘intelligent transmitters’’ which are able to exchange information, and to perform internal data processing and advanced functionalities. To make a contribution to solving these challenges, the mathematical framework of dynamic reliability is extended to handle the data and information which are processed and exchanged between systems components. Stochastic deviations that may affect system properties are also introduced to enhance the modelling of failures. A formalized Petri net approach is then presented to perform the corresponding reliability analyses using numerical methods. Following this formalism, a versatile model for the dynamic reliability modelling of digital-based transmitters is proposed. Finally the framework’s flexibility and effectiveness is demonstrated on a substantial case study involving a simplified model of a nuclear fast reactor. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Dynamic reliability Probabilistic dynamics Digital-based system Intelligent transmitter Smart sensor Petri net Piecewise-deterministic Markov process

1. Introduction Probabilistic risk assessments (PRA) (or ‘‘probabilistic safety assessments (PSA)’’) provide a general framework for managing risks linked to engineering systems, notably in nuclear power plants, chemical industries, and aerospace. The main purpose of PRA is to identify the possible accidental scenarios (i.e. sequences of events that could lead to an accident or undesired event), to rate their consequences (i.e. severity), and to assess their likelihoods (i.e. probabilities of occurrence, or ‘‘frequencies’’). An essential challenge is to deal with the system complexity, that is, to treat system interactions both efficiently and realistically. In the 1970s, the event tree/fault tree methodology was introduced for PRA [1], focusing on system components and their relationships. The event trees are inductive techniques which develop the possible events following an initiating event, taking the fulfilment or unfulfillment of safety functions (operations and

n Corresponding author at: Universite´ de Technologie de Troyes (UTT), Institut Charles Delaunay (ICD) & UMR CNRS 6279 STMR, 12 rue Marie Curie, BP 2060, 10010 Troyes Cedex, France. Tel.: + 33 3 25 75 96 90; fax: + 33 3 25 71 56 99. E-mail address: fl[email protected] (F. Brissaud).

0951-8320/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ress.2010.12.014

failures of systems, including safety barriers) into account. To complete the latter, fault trees are deductive techniques which illustrate and express safety functions (depicted as top events) as combinations of basic events, using logic gates. To perform analyses, the mathematical framework is then based on Boolean algebra. The event tree/fault tree methodology is now a widely used and accepted tool for risk management. In particular, it has proved its ability to model functional dependencies between system components (for example, through inclusion of common cause failures) [2]. However, these approaches have not been designed to deal with changes in event relationships (in sequences or in combinations) evolving according to time, for example due to interactions with process, environment, and humans [3]. They are therefore described as ‘‘static.’’ Other ‘‘static’’ models commonly used in reliability analyses include, for example, the reliability block diagrams. In the late 1980s, early 1990s, dynamic reliability (or ‘‘probabilistic dynamics’’) methods were developed to explicitly handle the influence of time, process dynamics, and human actions, on system operations and failures, and accidental scenarios. The first dynamic approach was denoted Dynamic logical analytical methodology (DYLAM) [4,5] and uses time discretization to simulate all the possible event sequences according to the evolution of

794

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Nomenclature set of non-negative integers set of real numbers t0 initial time, with t0 A R and t0 Z0 t time variable, with t A R and t Zt0 N number of components state variables, with N A N M number of process variables, with M A N L number of data variables, with L A N Q number of deviation variables, with Q A N i(t) vector of components state variables at time t, with iðtÞ A NN x(t) vector of process variables at time t, with xðtÞ A RM y(t) vector of data variables at time t, with yðtÞ A RL – y(t) vector of data variables at time t–e, that is –y(t)¼y(t–e), when e tends to 0 + (i.e. e tends to 0 with e 40), with yðtÞ A RL e(t) vector of deviation variables at time t, with eðtÞ A RQ p(ik-il9x(t), y(t), e(t), t) components transition rate from state ik to state il at time t, given x(t), y(t), and e(t) lik(x(t), y(t), e(t), t) total components transition rate from state ik at time t, given x(t), y(t), and e(t) Fik(t 9 x(t), y(t), e(t), t) probability of leaving components state ik in time interval [t, t+ t), given i(t) ¼ik, x(t), y(t), and e(t) [i] vector of all the possible values of vector i(t), for any time t, (i.e. each component of vector [i] is a distinctive vector of possible components state) Ei number of possible values of vector i(t), (i.e. dimension of vector [i]) P[i],i(t)(x(t), y(t), e(t), t)T  [i] product of vectors used to determine randomly i(t + Dt) according to transition rates, given i(t), x(t), y(t), and e(t) x0i(t)(x(t), e(t), t) expression of derivatives of process variables at time t, given i(t), x(t), and e(t), with xiðtÞ u ðxðtÞ, eðtÞ, tÞ A RM

N R



yi(t)(x(t), –y(t), e(t), t) expression of data variables at time t, given i(t), x(t), –y(t), and e(t), with yiðtÞ ðxðtÞ, yðtÞ, eðtÞ, tÞ A RL dEi(t)(x(t), e(t), t, t+ dt) vector of random increments of the deviation variables in time interval [t, t + dt], given i(t), x(t), and e(t), with dEiðtÞ ðxðtÞ, eðtÞ, t, t þ dtÞ A RQ xdEi(t)(de, x(t), e(t), t, t +dt) probability density function of dEi(t)(x(t), e(t), t, t + dt), using the variable of integration deA RQ , given i(t), x(t), and e(t) (i(t), x(t), y(t), e(t)) description of the complete system state at time t D safety domain of the complete system state, often defined by a boundary of process variables Gx such as the complete system is in a safe state at time t if and only if x(t) r Gx v(t) any vector at time t v(t)k kth component of vector v(t) vk specific value of vector v(t), indexed by k [v] vector of all the possible values of vector v(t), for any time t, (i.e. [v]k ¼vk for k¼1, 2, y, Ev, and vk avl for kal) Ev number of possible values of vector v(t), (i.e. dimension of vector [v]) [v]k kth component of vector [v] (with k¼1, y, Ev) [vk] vector of all the possible values of v(t)k, for any time t vk 7 ev interval defined by [vk–ev, vk + ev] v(t)T transpose of vector v(t) v(t)(S) subset (S) of components of vector v(t) v(t)(S),k kth component of vector v(t)(S) v(t)\(S) vector v(t) without subset (S) of its components Dt time step I1(A) indicator function defined by I1(A)¼1 if assertion A is true, and 0 otherwise Pr[E] probability of event E Pr[E9G] probability of event E, given event G



tures [16–18], and equivalent dynamic Bayesian networks [19].

applied to the well known ‘‘holdup tank problem’’ [25]: DDET [3], GO-FLOW methodology [3], Petri nets [22,24], Markov models [25,26], PDMP [28], and Monte Carlo methods [31]; or to systems taking part of nuclear reactors: DDET [4–9], ESD [13], dynamic fault trees [16], DFM [36–38], Markov models [36–38], and Monte Carlo methods [30,32]. The formal mathematical framework of the dynamic reliability was established under the name of Continuous event tree (CET) theory [39], introducing two sets of variables used to define the complete system state:

 The Dynamic flowgraph methodology (DFM) [20,21].  Petri nets with stochastic [22,23] and coloured characteristics

 State (discrete) variables of the system components (e.g.

process variables. Then, several other methods have emerged for dynamic reliability modelling:

 Discrete dynamic event trees (DDET), including DYLAM [4–6], DETAM [7], ADS [8], MCDET [9], and ADAPT [10].

 Event sequence diagrams (ESD) [11], with an extended framework for dynamic reliability [12,13].

 The GO-FLOW methodology [14,15].  Fault trees and reliability block diagrams with dynamic fea-



 

[24]. Markov models, combined with the cell-to-cell mapping technique (CCMT) to discretize time and other variables [25,26], or with the extension denoted (time-)continuous cell-to-cell mapping technique (CCCMT) [27]. Models which are directly based on Piecewise-deterministic Markov processes (PDMP), for performing exact Monte Carlo sampling methods [28], and for using deterministic algorithms [29]. Direct Monte Carlo approaches [30–32], discrete event simulations [3], and simulations supported by a stochastic hybrid automation [33].

Most of these methods have been presented, compared and discussed by Aldemir et al. [34], Siu [3], and Labeau et al. [35], and

operational states, failures).

 Physical (continuous) variables of the process (e.g. level, pressure, temperature). Changes in the components state variables were assumed stochastic, defined by transition rates, and dependent on process variables. On the other hand, the evolution of the process variables is determined by a set of first order (non-stochastic) differential equations, according to the components states. Components state and process variables are therefore mutually dependent, and follow a piecewise-deterministic (Markov) process (PDP) [40,41]. Under Markovian assumptions, Chapman–Kolmogorov equations have then been used to express the probability density function of the complete system state according to time [39,42]. Several works

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

were done to handle this probability distribution, by using method of moments [43], or by characterizing marginal distributions [44,45], and numerical solutions were provided [46]. The mathematical framework was also extended to introduce human operations by additional variables [47]. More recently the use of digital-based systems for safetyrelated applications has also introduced new issues for dynamic reliability modelling, especially due to interactions between system components [48,36,38], and when using communication networks [49–51]. Since the late 1990s, these considerations are increasing, notably in nuclear power plants where the analogue instrumentation is becoming obsolete and is being replaced by digital systems [48,52,36,38]. (It should be noted that similar questions have been raised early on in the aerospace industry for their large body of software-based control systems [53].) In addition to practical advantages and cost reductions, digital systems also offer the potential to improve safety and reliability by the use of their advanced features [54,48,55]. Typical examples of digital-based systems include the ‘‘intelligent transmitters’’ [56] of which features are shown in Table 1. In particular, these transmitters may cooperate by processing and exchanging data and information in order to improve availability and/or safety [65,66]. The dynamic reliability modelling of digital-based transmitters, and of systems integrating digital-based transmitters, is the key topic of the present paper. The objectives are to propose a framework for handling the dynamic interactions between system components, including ‘‘intelligent transmitters,’’ and between these components and process dynamics, and for performing safety assessments. Foremost, note that although there is a broad consensus on the need for dynamic reliability methods [67], and especially for complex and digital-based systems, these approaches have not yet penetrated the arena of industrial applications [35]. Given reasons are the theoretical flavour and the lack of a generic platform for performing these analyses. In particular, almost all the corresponding software tools are totally application specific [35]. In addition, since a general theory has to be observed in order to compare the different approaches [42], it has been considered that the hereafter proposed developments in dynamic reliability should:

 Be compatible with the mathematical framework of the Continuous event trees defined by Devooght and Smidts [39].

 Provide a generic approach, inherently appropriate to a wider field of applications with minimal relevant adjustments.

To this end, the present paper first extends the mathematical framework of dynamic reliability in Section 2, in order to

795

deal with the specific characteristics of digital-based systems. A formalized Petri net approach is then presented in Section 3, which provides a generic framework to express the model of the system and to perform the related reliability analyses. Following this formalism, a versatile model for the dynamic reliability modelling of digital-based transmitters is also proposed. A comprehensive case study involving a nuclear fast reactor is then presented in Section 4, and dynamic reliability modelling and analyses are performed in Section 5, to illustrate the application of the proposed approach.

2. Mathematical framework for dynamic reliability 2.1. Problem formulation In addition to time t, four types of variables are used to describe the complete system state. The components state and process variables remain the same as those defined by Devooght and Smidts [39]. Data variables are introduced in order to deal with the characteristics of digital-based systems; and deviation variables extend the possibilities of modelling failures. All these variables are time-dependent. Moreover, components state variables are discrete and therefore depicted by a vector of integers, denoted i(t); and the process, data, and deviation variables are continuous and therefore depicted by vectors of reals, denoted x(t), y(t), and e(t), respectively; making a hybrid (discrete/continuous) system. The notations are given in the Nomenclature. The components state variables i(t) represent the structure (configuration) of the system according to the physical states (operational, degraded, and failed) of its components and of human operations (e.g. opening or closing a valve). The state of any system component (e.g. operating, degraded, and failure modes) is described by an integer n, or by a set of integers (Sn), which are arranged in vector i(t), with iðtÞ A NN . (For example, a system that is composed of N components, each of these component described by state ni such as ni ¼ 1 if component i is operating and ni ¼0 otherwise, with i¼1, y, N, may be described by vector i(t)¼(n1, y, nN)T.) The components state variables may evolve both deterministically and stochastically, depending on process and data variables (e.g. an actuator is controlled by a signal (i.e. a data variable), and has a failure rate that depends on the temperature (i.e. a process variable)), and deviations (e.g. after a certain level of degradation (i.e. a deviation variable), a component’s transition occurs from a degraded mode into a fully nonoperating mode). The components transition rate from state ik to state il at time t, given process, data, and deviation variables (i.e. x(t), y(t), and e(t)), is denoted p(ik-il9x(t), y(t), e(t), t).

Table 1 What is an ‘‘intelligent transmitter’’? Feature

Description

Sensor or transmitter?

Sensor systems that combine data acquisition from physical or chemical properties with internal data processing in order to transmit an elaborate signal are appropriately referred to as ‘‘transmitters’’ instead of ‘‘sensors’’ (see e.g. [57])

Smart or intelligent?

Such systems are often described as ‘‘smart’’ if they incorporate signal conditioning and processing functions, carried out by embedded microprocessors [58–60], and ‘‘intelligent’’ depending on further functionalities involved in the external system functions [61,59], but this is not an unbreakable rule [56]

What functionalities?

Error measurement correction [58,60], self-adjustment [62], self-diagnosis and validation [63], online reconfiguration [64], and digital bidirectional communication [61,57]

What benefits?

More accurate measurement results [60], cost reductions [58], and use facilities

Dependable?

‘‘Intelligent transmitter’’ features have both pros and cons for reliability, maintainability, and safety [55,56]

796

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

These transition rates may take human operations into account (e.g. when a certain pressure level is detected through a signal (i.e. a data variable), an operator opens a valve (i.e. a component’s state is changed)). Then, the total components transition rate from state ik is X k lik ðxðtÞ, yðtÞ, eðtÞ, tÞ ¼ pði -il 9xðtÞ, yðtÞ, eðtÞ, tÞ ð1Þ il a ik

The transitions between components states are assumed instantaneous. When the components state at time t is i(t)¼ ik, the probability that the components leave this state before time t + t is therefore Fik ðt9xðtÞ, yðtÞ, eðtÞ, tÞ  Z t  lik ðxðt þuÞ, yðt þ uÞ, eðt þuÞ, t þuÞ  du ¼ 1exp 

ð2Þ

0

The process variables x(t), (xðtÞ A RM with M the number of process variables), represent the physical variables involved in the system dynamics (e.g. pressure, temperature, level, and volume). The process variables evolve deterministically, given components states, and with deviation variables as parameters (e.g. a level in a tank is determined by the configuration of the valves (i.e. a components state), and an amount of leakage (i.e. a deviation variable)). The evolution of process variables may then usually be defined by a set of first order (non-stochastic) differential equations: d xðtÞ ¼ xiðtÞ u ðxðtÞ, eðtÞ, tÞ dt

ð3Þ

The data variables y(t), ðyðtÞ A RL with L the number of data variables), represent any data or information which are computed, stored, and/or exchanged between system components (e.g. commands, measurement results, and diagnostic information). By nature, the data variables do not directly affect process and deviation variables, but may be used to change components state (e.g. a signal commands the activation of a unit). The data variables may usually be expressed as a function of process and deviation variables, given components state (e.g. when a transmitter is in a (non-fully) operating mode of operation (i.e. a component’s state), its measurement result depends on the quantity to be measured (i.e. a process variable), and drifts (i.e. deviation variables)) and may be initiated by human operations (e.g. a command is transmitted by a user interface). They may also depend on previous values of data variables (e.g. a result is computed using stored data; a signal is locked to a current value) which are denoted –y(t), defined by –y(t)¼ y(t  e) with e which tends to 0 + (i.e. e tends to 0 with e 40). Because the time t is given as an explicit parameter, –y(t) is able to include values of data variables at any time up to time t, using the relevant number L of data variables (e.g. a function is used to change components of y(t) only if time t fulfils specific conditions; an example is given in Appendix A). Finally the data variables are expressed as follows:



yðtÞ ¼ yiðtÞ ðxðtÞ, yðtÞ, eðtÞ, tÞ

ð4Þ

The deviation variables e(t), (eðtÞ A RQ with Q the number of deviation variables), represent continuous deviations or errors in system properties (e.g. system degradations and drifts). The deviation variables evolve stochastically, given components state, and with process variables as parameters (e.g. when a valve is closed (i.e. a component’s state), a leak follows a random distribution influenced by the pressure (i.e. a process variable)). A stochastic process [68] may then usually represent the evolution of deviation variables, defined by random increments in time

interval [t, t +dt]: eðt þ dtÞ ¼ eðtÞ þ dEiðtÞ ðxðtÞ, eðtÞ, t, t þ dtÞ

ð5Þ

with the random increments dEi(t)(x(t), e(t), t, t + dt) which follow a distribution with probability density function xdEiðtÞ ðde, xðtÞ, eðtÞ, t, t þ dtÞ (de is the variable of integration, used to represent the realisations of the random variable dEi(t)(x(t), e(t), t, t+ dt)). In Eqs. (3)–(5), which define the continuous variables, the components state (i.e. i(t)) is given as an index. In fact, one set of equations has to be theoretically defined for each state of components (i.e. system configuration). In practice, indicator functions such as I1(A), defined by I1(A)¼1 if assertion A is true, and 0 otherwise, may be however used to have only one set of equations for the continuous variables. Assertion A may then refer to the state of all the system components or, in most cases, to only a part of system components, given in a subset (S) of components of vector i(t), denoted i(t)(S). In Eqs. (1)–(5) it is also possible to explicitly specify a vector of parameters, usually denoted a [42,35], which is, for example, useful when performing uncertainty analyses [69,70]. The complete state of the system is then described by (i(t), x(t), y(t), e(t)). To perform risk assessment, a safety domain D is therefore defined such as the complete system is in a safe state if and only if (i(t), x(t), y(t), e(t))AD. Usually D only refers to process variables, that is, a hazardous event is assumed if and only if some process variables (e.g. temperature, pressure) reach certain values. D is therefore defined by a boundary of process variables Gx such as the complete system is in a safe state at time t if and only if x(t)r Gx. Finally reliability analyses aim to assess the probability, according to time, that the complete system state leaves this safety domain, and more generally, risk assessment consists in investigating the corresponding scenarios. 2.2. Numerical solution Even without data and deviation variables, an analytical solution for the probability of the complete system state at time t has appeared unmanageable, apart from very simple cases [44], and only numerical solutions have been provided [44,46]. Moreover, the use of numerical approaches circumvents the need to discretize the continuous variables a priori, as required in several dynamic reliability methods, including DFM and Markov models with CCMT. In fact, the discretization of continuous variables such as process variables may lead to drawbacks in terms of modelling size capability versus analysis accuracy [35]. These considerations further increase when data and deviation variables are added and, for the latter variables, it is also not always obvious to know how to a priori define efficiently ranges of evolution. Discussions about numerical methods which can be applied to these topics can be found in the literature (e.g. [71–73]). Since transition rates and differential equations are used to describe the system evolution (cf. Eqs. (1)–(5)), the numerical method adopted in the present paper is based on Taylor series and finite differences. A time step Dt is then introduced. It should be small enough to assume that variables i(t), x(t), y(t), and e(t) are constant in any time interval [t, t + Dt), without loss of accuracy. These variables at time t + Dt (i.e. i(t + Dt), x(t+ Dt), y(t + Dt), and e(t + Dt)) can then be determined according to their values at time t (i.e. i(t), x(t), y(t), and e(t)), using Eqs. (1)–(5). In particular, a components state transition that occurs between time t and time t+ Dt is considered to occur exactly at time t+ Dt. In the same way, the evolutions of the process, data, and deviation variables between time t and time t + Dt are considered to occur as ‘‘jumps’’ exactly at time t + Dt. It is then possible to approximate the probability that the components remain in their current state at time t, denoted

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

i(t)¼ik, up to time t + Dt, that is i(t+ Dt)¼ik, according to Eq. (2), and using the Taylor series (here, a first order is used, other orders can also be expressed if Dt is not small enough to yield accurate approximations): Pr½iðt þ DtÞ ¼ ik 9iðtÞ ¼ ik , xðtÞ, yðtÞ, eðtÞ, t ¼ 1Fik ðDt9xðtÞ, yðtÞ, eðtÞ, tÞ  1-Dt  lik ðxðtÞ, yðtÞ, eðtÞ, tÞ

ð6Þ

Similarly the probability that the components leave their current state at time t, denoted i(t) ¼ik, for another specific state at time t + Dt, denoted i(t + Dt)¼il, with ik ail, can be approximated by l

k

k

Pr½iðt þ DtÞ ¼ i ai 9iðtÞ ¼ i , xðtÞ, yðtÞ, eðtÞ, t

 Dt  pðik -il 9xðtÞ, yðtÞ, eðtÞ, tÞ

ð7Þ

Note that deterministic (certain) transitions can therefore be modelled using transition rates that are equal to 1/Dt. To determine the components state at time t + Dt (i.e. i(t+ Dt)) according to Eqs. (6) and (7), the following procedure can be used: i. First, vectors are defined to represent the possible values of components states and the corresponding probabilities of transitions: let [i] be a vector defined by all the possible values of vector i(t) (for any time t), that is, [i]k ¼ ik for k¼1, 2, y, Ei, ([i]k is the kth component of vector [i]), with Ei the total number of possible values of vector i(t) (i.e. the total number of possible components states), and ik ail for kal; given the system variables at time t (i.e. i(t), x(t), y(t), and e(t)), let q[i],i(t)(x(t), y(t), e(t), t) be a vector defined by the probabilities of transitions from state i(t), that is q[i],i(t)(x(t), y(t), e(t), t)k ¼ Dt  p(i(t)-ik9x(t), y(t), e(t), t) if i(t)aik, and q[i],i(t)(x(t), y(t), e(t), t)k ¼1–Dt  lik(x(t), y(t), e(t), t) if i(t) ¼ik, for k¼1, 2, y, Ei, (q[i],i(t)(x(t), y(t), e(t), t)k is the kth component of vector q[i],i(t)(x(t), y(t), e(t), t)). ii. Then, the simulation of components state at time t + Dt is obtained through the use of two random variables: let R[i],i(t)(x(t), y(t), e(t), t) be a discrete random variable with sample space {1, 2, y, Ei} and probability mass function defined by q[i],i(t)(x(t), y(t), e(t), t), that is, Pr[R[i],i(t)(x(t), y(t), e(t), t)¼k]¼ q[i],i(t)(x(t), y(t), e(t), t)k for k¼1, 2, y, Ei; and let P[i],i(t)(x(t), y(t), e(t), t) be a vector with all its components equal to 0, except for its component indexed by the realisation of R[i],i(t)(x(t), y(t), e(t), t) that is equal to 1. iii. Finally iðt þ DtÞ  P [i],iðtÞ ðxðtÞ, yðtÞ, eðtÞ, tÞT U½i

ð8Þ

797

Eqs. (6)–(11) show that the complete state of the system at time t + Dt, that is (i(t+ Dt), x(t + Dt), y(t + Dt), e(t + Dt)), can be fully determined according to the complete state of the system at time t, that is (i(t), x(t), y(t), e(t)), according to deterministic and stochastic evolutions. The system therefore follows a piecewisedeterministic (Markov) process (PDP) [40,41]. The following developments aim at expressing the probability of the complete system state at time t+ Dt, according to the complete system state at time t. Because Eqs. (9)–(11) provide numerical approximations of the results, the following constant vectors are introduced: ex A RM , ey A RL , and ee A RQ , to express the complete system state at time t by (i(t), x(t), y(t), e(t))¼ (ik, xk 7 ex, yk 7 ey, ek 7 ee), with vk 7 ev which represents the interval defined by [vk–ev, vk + ev]. Moreover, the vectors ex, ey, and ee, should be small enough to assume the following approximation: Pr½ðiðt þ DtÞ, xðt þ DtÞ, yðt þ DtÞ, eðt þ DtÞÞ ¼ ðik , xk 7 ex , yk 7 ey , ek 7 ee Þ9ðiðtÞ, xðtÞ, yðtÞ, eðtÞÞ ¼ ðil , xl 7 ex , yl 7 ey , el 7 ee Þ  Pr½ðiðt þ DtÞ, xðt þ DtÞ, yðt þ DtÞ, eðt þ DtÞÞ ¼ ðik , xk 7 ex , yk 7 ey , ek 7 ee Þ9ðiðtÞ, xðtÞ, yðtÞ, eðtÞÞ ¼ ðil , xl , yl , el Þ

ð12Þ

for any complete system states (ik, xk, yk, ek) and (il, xl, yl, el). The conditional probability that the complete system is in a specific state at time t+ Dt, denoted (i(t + Dt), x(t + Dt), y(t + Dt), e(t + Dt))¼(ik, xk 7 ex, yk 7 ey, ek 7 ee), given it was in a specific state at time t, denoted (i(t), x(t), y(t), e(t)) ¼(il, xl, yl, el), can then be approximated as follows: Pr½ðiðt þ DtÞ, xðt þ DtÞ, yðt þ DtÞ, eðt þ DtÞÞ ¼ ðik , xk 7 ex , yk 7 ey , ek 7 ee Þ9ðiðtÞ, xðtÞ, yðtÞ, eðtÞÞ ¼ ðil , xl , yl , el Þ  ½I1 ðil ¼ ik Þ  ð1Dt  lil ðxl , yl , el , tÞÞ þ I1 ðil a ik Þ  Dt  pðil -ik 9xl , yl , el , tÞ  I1 ðxk ex rxl þ DtUxiul ðxl , el , tÞ rxk þ ex Þ  I1 ðyk ey ryik ðxk , yl , ek , t þ DtÞ r yk þ ey Þ  ð2  99ee 99  xdE l ðek , xl , el , t, t þ DtÞÞ i

ð13Þ

with I1(A) the indicator function defined in the Nomenclature, and 99 ee 99 the Euclidian norm of vector ee. And the unconditional probability that the complete system is in a specific state at time t + Dt, denoted (i(t + Dt), x(t + Dt), y(t + Dt), e(t + Dt))¼(ik, xk 7 ex, yk 7 ey, ek 7 ee) is Pr½ðiðt þ DtÞ, xðt þ DtÞ, yðt þ DtÞ, eðt þ DtÞÞ ¼ ðik , xk 7 ex , yk 7 ey , ek 7 ee Þ X Z þ1 Z þ1 Z þ1  Pr½ðiðt þ DtÞ, xðt þ DtÞ, yðt þ DtÞ, eðt þ DtÞÞ il

1

1

1

k

The values of process variables at time t + Dt (i.e. x(t + Dt)) can be approximated according to Eq. (3), and using finite differences (with dt¼ Dt): xðt þ DtÞ  xðtÞ þ Dt  xuiðtÞ ðxðtÞ, eðtÞ, tÞ

ð9Þ

The values of deviation variables at time t + Dt (i.e. e(t + Dt)) are generated as paths of the stochastic process model defined by Eq. (5), using dt ¼ Dt: eðt þ DtÞ ¼ eðtÞ þ dEiðtÞ ðxðtÞ, eðtÞ, t, t þ DtÞ

ð10Þ

Once the components state, process and deviation variables are defined at time t+ Dt, the data variables at time t + Dt (i.e. y(t + Dt)) can be also determined according to Eq. (4), with e ¼ Dt and –y(t + Dt)¼y(t): yðt þ DtÞ  yiðt þ DtÞ ðxðt þ DtÞ, yðtÞ, eðt þ DtÞ, t þ DtÞ

ð11Þ

¼ ði , xk 7 ex , yk 7 ey , ek 7 ee ÞjðiðtÞ, xðtÞ, yðtÞ, eðtÞÞ ¼ ðil , xl , yl , el Þ del  dyl  dxl

ð14Þ

Starting from initial conditions at time t0, denoted (i(t0), x(t0), y(t0), e(t0)), it is therefore possible to assess the probability of the complete system state at time t, using an iterative process based on Eqs. (6)–(14). The reliability analyses then focus on the following probability: Pr½ðiðtÞ, xðtÞ, yðtÞ, eðtÞÞ A D ¼ Pr½xðtÞ r Gx 

ð15Þ

It is obvious that the number of possible scenarios, according to all the variables, may be unmanageable. In the following, a Monte Carlo approach is therefore adopted to compute Eq. (15) by iteration of Eqs. (6)–(11). Unlike other commonly used Monte Carlo methods, the proposed approach does not simulate the set points (i.e. next time instants where variables are changed), but simulates the

798

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

evolution of the variables in every time step Dt (the set points are therefore fixed and equal to t0 +k  Dt with k¼1, 2, y). The evolution of the system variables forms a scenario which can be classified according to its outputs. By repeating these random samplings, the probability at time t, defined by Eq. (15), may be assessed by the relative frequency of occurrence of event x(t)r Gx. Direct simulations in programming languages could therefore be used to perform these reliability analyses. However, such approaches usually suffer from a lack of clarity and transparency, and changes in the modelling are therefore restricted to the programmer. To have a visual tool based on a common approach well known by safety and reliability engineers, a Petri net formalism is then proposed in the next section, providing a generic framework to perform the dynamic reliability analyses using the above-presented numerical methods.

hybrid Petri nets [82] and differential Petri nets [83]. An additional advantage of Petri nets is the number of software tools available [84], including high quality free tools. It is not in the scope of the present paper to propose a new Petri net extension. Stochastic and coloured extensions are assumed to provide the two general properties useful for dynamic reliability. However, the specific manner in which these concepts are used needs to be defined to maximise the benefits of modelling and analysis. In the present paper, a Petri net formalism is therefore presented, using stochastic and coloured properties, to provide a generic framework:

 To flexibly model the dynamic reliability of a system with the help of a visual and easy-to-handle interface.

 To simulate the evolution of the complete system state and perform reliability analyses, using numerical methods and Monte Carlo simulations.

3. Formalism for dynamic reliability modelling and analyses 3.1. Petri net formalism On the one hand, dynamic reliability has to deal with transitions based on set points (time and variable conditions), competing events, and simultaneous transitions and sequences [35]. On the other hand, Petri nets provide an interesting tool for describing and studying systems that are characterised as being concurrent, asynchronous, distributed, parallel, nondeterministic, and/or stochastic [74]. The use of Petri nets for dynamic reliability is therefore natural, and especially when dealing with digitalbased systems. Moreover, the stochastic and coloured Petri nets [75,76] are promising extensions for modelling dynamic systems [77], notably for risk analysis [78]. For example, stochastic [22], and coloured [24] Petri nets have been already applied to the ‘‘holdup tank problem,’’ and for modelling cooperating transmitters [65,66]. The theory of Petri nets is well developed and continuously updated, as shown by the numerous new named extensions of Petri nets that regularly appear in the scientific literature. Several of these extensions offer interesting properties to deal with interactions between discrete and continuous variables: hybrid Petri nets [79] where both real and integer numbers can be used to mark places; extended coloured Petri nets [80] which model ‘‘dynamic colours’’; and dynamically coloured Petri nets [81] which explicitly link Petri nets and PDP. The computation of differential equations may also take advantage of high-level

In the proposed approach, each place of the Petri net is associated to one set of variables, and vice-versa. The number of places then increases linearly with the number of variables, which avoids any combinatorial explosion in the modelling. According to the nature of the variables (continuous or discrete, stochastic or deterministic), different representations are used for graphical convenience, as shown in Fig. 1. The values of the variables are given by the (coloured) tokens, with real or integer numbers (according to the places) inside the corresponding places, and are changed by the transitions. Each place therefore always contains one and only one token, and thus, all the transitions are always enabled. Guards are then used for each transition, and denoted sj[Dt], which means that the transition is fired at each time instant sj + k  Dt, with k¼0, 1, 2, y Moreover, the values specified on the arcs (when relevant) do not refer to a number of tokens, but to the values of variables which are handled. Each transition of the Petri net is associated to one specific place, denoted ‘‘managed place.’’ This place is linked to the transition by an input arc (from the place to the transition), which means that the variables depicted in this place are changed by the transition (‘‘the token is removed from the place’’); and linked to the same transition by an output arc (from the transition to the place), which attributes to the variables their new values according to an expression (‘‘a new token is deposited in the place’’). The latter may depend on the previous values of the variables (handled by the input arc) and also on variables from other places. These latter places are denoted ‘‘dependence places’’

Fig. 1. Petri net elements.

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

and are linked to the corresponding transition by bi-directed arcs, which means that the values of their variables may be used by the transition, but are not changed. The input arcs (from the ‘‘managed places’’ to their associated transitions), and the bi-directed arcs (between the transitions and their ‘‘dependence places’’), therefore do not have to be valued because they always provide the current values of the variables depicted in the linked place. On the other hand, the output arcs (from the transitions to their ‘‘managed places’’) have to specify the values that are given to the variables of the ‘‘managed places’’ after each solicitation. Depending upon the ‘‘managed places’’ considered, these values may be random variables. Contrarily to the ‘‘classical’’ stochastic Petri nets, the stochastic aspects therefore do not relate to time instants, but to values of variables (i.e. the token ‘‘colour’’). The time, as a variable, is modelled by a place. Besides, all the transitions are fired at deterministic time instants (specified by the transition guards). These Petri nets can therefore be classified as an ‘‘untimed stochastic (and coloured) Petri net.’’ A generic Petri net for the dynamic reliability modelling of systems is shown in Fig. 2, using the elements described in Fig. 1. The five types of variables defined in Section 2.1 are modelled, and depicted using different varieties of places. In Fig. 2, the variables are represented by vectors (except for time t which is a scalar). For more detailed models, and to treat dependencies more efficiently, it is also possible to split these vectors into subsets with one place for each (subset), and to handle them separately (cf. the example of digitalbased transmitters presented in Section 3.2). Each transition is fired at every time step Dt, changing the variables modelled by the corresponding ‘‘managed place,’’ according to the equations given in Section 2.2 and specified on the output arcs. For example, the transition that changes the time variable t is not linked to any ‘‘dependence place’’ and is simply used to increment time t by Dt at each solicitation. On the other hand, the components state variables i(t) are changed by a random variable which depends on all the other variables. In each time interval [t, t+ Dt), all the values of variables at time t + Dt are computed following a specific order defined by the sj (cf. top of Fig. 2), and using the values of variables at time t, according to the equations given in Section 2.2. In particular, x(t + Dt) and e(t + Dt) both depend on x(t) and e(t). To avoid

799

‘‘loosing’’ the values of x(t) (resp. e(t)) after computing x(t+ Dt) (resp. e(t + Dt)), and still be able to compute e(t + Dt) (resp. x(t + Dt)), ‘‘meta-transitions’’ are introduced and shown in Fig. 3. They are used to first compute the evolutions of variables in time interval [t, t + Dt] (i.e. x0 i(t)(x(t), e(t), t) and dEi(t)(x(t), e(t), t, t + Dt)), storing them as additional variables, and then to change the system variables. A ‘‘meta-transition’’ has therefore a double guard denoted (sj, sk)[Dt], which means that the evolution of the variables is computed at each time instant sj + k  Dt, and the variables of the ‘‘managed place’’ at each time instant sk + k  Dt, with k¼0, 1, 2, y and sj osk (cf. Figs. 2 and 3). To guarantee that Eqs. (6)–(11) are well computed, a specific order for the sj is required (cf. top of Fig. 2): the evolutions of process and deviation variables are computed first (s1 and s2) because they do not change variables that are required by other equations; then the

Fig. 3. ‘‘Meta-transition’’ for Petri net.

Fig. 2. Generic Petri net for the dynamic reliability modelling of systems.

800

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Fig. 4. Model for dynamic reliability of digital-based transmitters.

components states variables (s3) because they require all the variables at the current time; the time t variable (s4); the process and deviation variables (s5 and s6) because they now depend only on their current values and evolutions previously computed (cf. Eqs. (9) and (10)); and finally the data variables (s7) because they require all the other variables at the next time instant (cf. Eq. (11)). To deal efficiently with the dependencies between subsets of variables, such specific orders may also be used between variables of the same type, when handling them separately. In addition, the ‘‘meta-transitions’’ make available additional variables (the evolutions of process and data variables) which can, for example, improve the modelling when the evolution of variables depends on derivatives of other variables.

3.2. Model for digital-based transmitters The ‘‘intelligent transmitters’’ take advantage of digital technologies to perform internal data processing and functionalities such as error measurement correction, self-adjustment, self-diagnoses and validation, online reconfiguration, and digital bi-directional communication [56]. In particular, they are able to exchange information and data (e.g. measurement results and diagnostic information) in order to improve their operations [65]. A versatile model is proposed in Fig. 4 for the dynamic reliability modelling of such digital-based transmitters, based on the framework given in Section 2, and the formalism presented in Section 3.1. This model does not pretend to be exhaustive. It is founded on previous published research on

‘‘intelligent transmitters’’ [56], and on the features given in Table 1. It should then be adapted for the proper systems modelling. If need be this modelling scheme can and should be revised as and if further information becomes available. To make the model clearer, the place that represents the time t variable has been omitted because it can potentially be linked by bi-directed arcs to almost all of the places. Moreover, the specifications of the output arcs and the transition guards are not given either, because these considerations are out of the scope of such a model. On the other hand, the focus is placed on the data variables used to represent information and data which are handled by the transmitter, and exchanged between the transmitter and other system components (including other transmitters). Subsets of data variables from vector y(t) are specified by indices in parentheses (cf. Fig. 4). For example, y(t)(T) contains all data variables that may be modified by the transmitter, and is itself comprised of subsets explained hereafter. On the contrary, y(t)\(T) contains all other data variables, which include the data and information received from other system components and transmitters, and human commands. Vector x(t) includes the quantities to be monitored by the transmitter, i(t)(T) gives the functional or dysfunctional state (operating, degraded, and failure modes, ‘‘detected’’ or not by self-diagnoses) of the transmitter, and e(t)(T) represents transmitter deviations such as drifts. As ‘‘central information,’’ the process data of the transmitter are given by y(t)(P) (note that the ‘‘process data’’ are part of data variables and should not be confused with process variables). They represent all the parameters that are computed by the

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

transmitter, in order to provide any information about its functional or dysfunctional state. These process data may be directly obtained through the transmitter state (i.e. i(t)(T)), for example if the transmitter is in a ‘‘detected’’ failure mode, and also using data from other system components (i.e. y(t)\(T)), such as measurement results and diagnostic information from other transmitters. A feedback on the measurement results of the transmitter (i.e. y(t)(M)) may also be used, for example to compare or compute these results with the measurement results from other transmitters. Finally relevant data may be stored (i.e. y(t)(D)), for example to assess the coherence of current information with previous results. Using the process data, if the transmitter state is assumed ‘‘inappropriate,’’ transmitter reconfigurations may be commissioned (i.e. y(t)(R)), changing the transmitter state (i.e. i(t)(T)), for example to activate a cold stand-by of the transmitter or to change a measure bandwidth. The measurement results of the transmitter (i.e. y(t)(M)) are then function of the transmitter state (i.e. i(t)(T)), the quantities to be measured (i.e. x(t)), and possible errors such as drifts (i.e. e(t)(T)). Moreover, parameters of correction (i.e. y(t)(C)), defined using the process data, may be used, for example to offset drifts. In addition to the measurement results, the transmitter is also able to transmit other data such as diagnostic information (i.e. y(t)(I)) to other system components, for example, a degree of confidence in the transmitter results. Depending on the specific applications, other relationships between elements of Fig. 4 (by the use of bi-directed arcs) may be added. For example, the stored data, the measurement results, and the data from other system components (i.e. y(t)(D), y(t)(M), and y(t)\(T)) may be handled only through the process data (assumption made in Fig. 4), or may also be directly used for the computation of other variables (such examples are used in the case study developed in Section 4).

4. Case study: description of a system with digital-based transmitters 4.1. The Europa fast reactor The case study described in this section is the primary circuit of the Europa fast reactor. This system has been proposed as a benchmark problem on accident sequences [85]. Several dynamic reliability analyses have been also performed on this application, using DYLAM [4], Monte Carlo methods [30], and the ESD approach [13]. A comprehensive description of the original system is provided by Smidts and Devooght [30]. In the present paper, the physical variables have been simplified (the effect of the reactivity in normal operations has been neglected, in accordance with previous works which have shown the incidental impact of the reactivity on safety functions [13]) to allow us to more specifically focus on other aspects. In particular, new transmitter features have been introduced (data processing, bi-directional communications, error measurement corrections), as well as deviation variables. Some notations and failure rates have been also changed from the original case study. The simplified model of the primary circuit of the Europa fast reactor is shown in Fig. 5. This system is comprised of two channels (C1 and C2) where sodium is introduced as a coolant by a pump (PM). The lack of coolant, for example in the case of a pump failure, increases the temperatures in the channels and may yield hazardous events. The safety function studied in the following consists in the detection of a high sodium temperature in the channels, or a low sodium flow rate, to activate an emergency system shutdown. The sodium temperature in each

801

Fig. 5. Primary circuit of the Europa fast reactor (simplified model).

channel (T1 and T2) is therefore monitored by a transmitter (resp. m ST1 and ST2), which sends its results (resp. Tm 1 and T2 ) to a common controller (CT). Similarly the sodium (linear momentum) flow rate (G) is monitored by a third transmitter (SG), which sends its result (Gm) to a dedicated controller (CG). The two primary controllers (CT and CG) send their binary signals (resp. yCT and yCG) to a central controller (SDL). If a high temperature threshold (Tmax) is detected in at least one channel (i.e. yCT ¼1), or if a low flow rate threshold (Gmin) is detected (i.e. yCG ¼1), then the central controller (SDL) sends a binary signal (ySDL) that should activate an emergency shutdown (SCM), (i.e. ySDL ¼ 1). This safety device, commonly named SCRAM, consists of the insertion of control rods (under gravity) into the reactor core that quickly stop the nuclear reaction by absorbing neutrons. The full description of the system variables (components state, process, data, and deviation variables) is provided in the following sections. Beforehand, note that the pump (PM) may have a full failure, which implies a nil output flow rate, or can be subject to torque deviations (denoted dM), which implies a decreasing output flow rate. Moreover, the temperature transmitters (ST1 and ST2) may be subject to drifts (resp. denoted dD1 and dD2), which m implies measurement results (resp. Tm 1 and T2 ) that are randomly biased. ‘‘Intelligent transmitter’’ features are then used to perform error measurement corrections. To this end, the three transmitters (ST1, ST2, and SG) are able to store, exchange, and process data in order to detect potential drifts, and to introduce parameters of corrections (denoted Tc1 and Tc2) in measurement results. 4.2. Components state variables Eight system components are assumed: two mechanical components which are the pump and the SCRAM (PM and SCM), three transmitters (SG, ST1, and ST2), and three controllers (CG, CT, and SDL). The state of each of these components is represented using a finite integer, as defined in Table 2. The vector of components state variables at time t is therefore i(t) ¼(SPM(t), SSG(t), SCG(t), SST1(t), SST2(t), SCT(t), SSDL(t), SSCM(t))T, and iðtÞ A NN with N ¼8. In the following, and in accordance with the formalism presented in Section 3.1, each component of vector i(t) is modelled separately. The mode of ‘‘normal’’ operation (i.e. full operating mode) of a component is represented by a state variable that is equal to 1 i.e. when SPM(t) ¼1, the pump provides its function according to the nominal performances; when SSCM(t)¼1, the SCRAM is able to be activated; when SSG(t)¼1, SST1(t) ¼1 and SST2(t) ¼1, the

802

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Table 2 Components state variables. Component

State variable

Value

Description

Pump (PM)

SPM(t)

¼1 ¼0 ¼3

Normal operation: the pump torque is nominal Full failure: the pump torque is nil Partial failure: the pump torque is subject to deviations

Flow rate transmitter (SG)

SSG(t)

¼1 ¼0 ¼2

Normal operation: the results are correct Dangerous failure: the results are locked to current values Safe failure: the results are locked at Gmin or below

Flow rate controller (CG)

SCG(t)

¼1 ¼0 ¼2

Normal operation: the signal is correct Dangerous failure: the signal is locked to 0 Safe failure: the signal is locked to 1

Temperature transmitter i (STi), with i¼1, 2

SSTi(t)

¼1 ¼0 ¼2 ¼3 ¼4

Normal operation: the results are correct Dangerous failure: the results are locked to current values Safe failure: the results are locked at Tmax or above Dangerous drifts: the results are subject to negative drifts Safe drifts: the results are subject to positive drifts

Temperature controller (CT)

SCT(t)

¼1 ¼0 ¼2

Normal operation: the signal is correct Dangerous failure: the signal is locked to 0 Safe failure: the signal is locked to 1

Central controller (SDL)

SSDL(t)

¼1 ¼0 ¼2

Normal operation: the signal is correct Dangerous failure: the signal is locked to 0 Safe failure: the signal is locked to 1

SCRAM (SCM)

SSCM(t)

¼1 ¼0 ¼5

Normal operation: the SCRAM can be activated Full failure: the SCRAM cannot be activated Activation mode: the SCRAM is activated

Table 3 Components transition rates. State variable

From state

To state

Transition rate

Valuea [s  1]

SPM(t)

1 or 3 1

0 3

lPM,0(dM(t)) lPM,3

¼1  10  3  exp(dM(t)  5  10  5) ¼1  10  2

SSG(t)

1 1

0 2

lSG,0 lSG,2

¼2  10  3 ¼2  10  4

SCG(t)

1

0 or 2

lCG,0/2

¼1  10  5

SST1(t)

1, 3 or 4 1, 3 or 4 1

0 2 3 or 4

lST1,0(SST2(t), T1(t)) lST1,2(SST2(t), T1(t)) lST1,3/4

¼4  10  4  (1+ I1(SST2(t) ¼ 0))  r(T1(t)) ¼4  10  5  (1+ I1(SST2(t) ¼ 2))  r(T1(t)) ¼1.5  10  2

SST2(t)

1, 3 or 4 1, 3 or 4 1

0 2 3 or 4

lST2,0(SST1(t), T2(t)) lST2,2(SST1(t), T2(t)) lST2,3/4

¼4  10  4  (1+ I1(SST1(t) ¼ 0))  r(T2(t)) ¼4  10  5  (1+ I1(SST1(t) ¼ 2))  r(T2(t)) ¼1.5  10  2

SCT(t)

1

0 or 2

lCT,0/2

¼1  10  5

SSDL(t)

1

0 or 2

lSDL,0/2

¼1  10  6

SSCM(t)

1

0

lSCM,0(T1(t), T2(t), t)

1

5

lSCM,5(ySDL(t))

¼3.125  10  8  (t  t0)  r((T1(t) + T2(t))/2) ¼I1(ySDL(t) ¼1)  (1/e)b

a With the indicator function I1(A) defined in the Nomenclature, and the following function used to model the effect of temperatures on transition rates: r(T) ¼6.17  10  2  exp(5.21  10  3  T). b With e which tends to 0 + such as lSCM,5 tends to infinity (i.e. deterministic transition) when ySDL(t) ¼ 1.

transmitters provide correct results (neither locked nor subject to drifts); and when SCG(t)¼1, SCT(t)¼ 1, and SSDL(t)¼1, the controllers provide correct signals. When a component state variable is equal to 0, the corresponding failure mode is ‘‘dangerous,’’ that is, it may yield the occurrence of hazardous events or disable the safety function: the mechanical components are not able to perform their function; the transmitter results are locked to current values (and thus, they are not able to detect new threshold values), and the controller signals are locked to 0 (i.e. signals that do not command the SCRAM activation). On the contrary,

when a component state variable is equal to 2, the corresponding failure mode is ‘‘safe,’’ that is, it may yield a spurious activation of the SCRAM (i.e. a SCRAM activation while the temperatures and the flow rate are under the threshold values): the transmitter results are locked at values which are equal or either below (for flow rate) and above (for temperatures) the threshold values; the controller signals are locked to 1 (i.e. signals that command the SCRAM activation). ‘‘Safe’’ failure modes are not assumed for the mechanical components. Other values of components state variables, which are 3 and 4, correspond to ‘‘degraded’’ modes of

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

operation (i.e. partial failure mode): the pump torque is subject to deviations; the results of the temperature transmitters are subject to positive or negative drifts (overestimated or underestimated results, respectively). ‘‘Degraded’’ modes of operation are not assumed for the other components. Finally the activation of the SCRAM corresponds to a specific mode, which is denoted by a SCRAM state variable that is equal to 5. The state variables of the mechanical components directly affect process variables as described in Section 4.3. On the other hand, the state variables of the controllers and transmitters directly determine data variables defined in Section 4.4. Finally the effects of the ‘‘degraded’’ modes of operation are modelled using deviations variables as described in Section 4.5. The transition rates between the possible states of each component are given in Table 3. Note that some transition rates depend directly on time t (a SCRAM transition rate fits a Weibull distribution with a shape parameter equal to 2.0, a scale parameter influenced by the sodium temperatures, and a location parameter equal to t0), process variables (temperature transmitters and SCRAM transition rates are function of the temperatures T1(t) and T2(t)), deviation variables (a pump transition rate is function of the pump torque deviation dM(t)), information variable (for the SCRAM activation according to the central controller signal ySDL(t)¼1), and also on the states of other components (the states of the temperature transmitters are stochastically dependent on each other due to common cause failures). Finally a deterministic (certain) transition is used for the SCRAM activation (from SSCM(t)¼1 to SSCM(t)¼5), when the SCRAM is in the mode of ‘‘normal’’ operation and the central controller sends a signal equal to 1 (i.e. ySDL(t)¼1). This transition is modelled using a rate that tends to infinity, i.e. equal to 1/e, with e which tends to 0 + . (For the numerical analyses using the formalism presented in Section 3.1, e ¼ Dt.)

803

4.3. Process variables The process dynamics are determined by thirteen process variables which are modelled separately and defined hereafter by first order differential equations. The vector of process variables at time t is x(t)¼(o(t), G(t), Tc,1(t), Tc,2(t), T1(t), T2(t), P(t), C1(t), C2(t), C3(t), C4(t), C5(t), C6(t))T, with xðtÞ A RM and M¼13. The parameters used in the definition of the process variables are given in Table 4, and the initial conditions at time t ¼t0, which correspond to steady states of process variables (i.e. x(t) is constant over time t, if the deviation variables are nil and the components state variables correspond to the modes of ‘‘normal’’ operation) are given in Table 5. The angular speed of the pump, denoted o(t) [rad s  1], is given by d C  I ðSPMðtÞ a 0ÞdM ðtÞ  I1 ðSPMðtÞ ¼ 3ÞK  oðtÞ oðtÞ ¼ M 1 dt I

ð16Þ

with the pump state variable SPM(t) defined in Section 4.2 (cf. Table 2), the pump torque deviation variable dM(t) defined in Section 4.5, and the indicator function I1(A) defined in the Nomenclature. Other parameters are defined in Table 4. The sodium linear momentum flow rate, denoted G(t) [kg m  2 s  1], is given by d u  GðtÞ2 GðtÞ ¼ þ Cnt2  ðoðtÞ2 oðt0 Þ2 Þu  Gðt0 Þ dt Gðt0 Þ

ð17Þ

The fuel temperature in channel i, denoted Tc,i(t) [K], with i¼1, 2, is given by d R  wi  PðtÞðTc,i ðtÞTe Þ T ðtÞ ¼ dt c,i ti ðTc,i ðtÞÞ

ð18Þ

Table 4 Parameters of process variables. Parameter

Value

Unit

Description

CM K I u Cnt2 R w1 w2 Te t1(T)

60,000 10 10 –5.00  10  2 5.56  10  6 9.521  10  7 0.41175 0.58825 653 1.1223+ 1.5215  10  3  T 1.0471  10  6  T2 + 2.7476  10  10  T3 1.5714+ 2.1303  10  3  T 1.4661  10  6  T2 + 3.8470  10  10  T3 0.518867 0.741238 1629  8.3290  10  1  T 4.40  10  1 8.2100  10  5 7.4480  10  4 6.6150  10  4 1.3277  10  3 6.1480  10  4 1.8940  10  4 3.98  10  2 1.29  10  2 3.11  10  2 1.34  10  1 3.31  10  1 1.26 3.21

Nm kg m2 s  1 kg m2 s1 kg m  2 K W1 – – K s

Nominal pump torque Constant which models pump friction phenomena Pump moment of inertia Constant for flow rate evolution Constant for flow rate evolution Thermal resistance Proportion of power generated by C1, with w1 + w2 ¼ 1 Proportion of power generated by C2, with w1 + w2 ¼ 1 Sodium reactor inlet temperature Time constant which includes the fuel specific heat in C1

s

Time constant which includes the fuel specific heat in C2

m2 m2 J kg  1 K  1 s1 – – – – – – s s1 s1 s1 s1 s1 s1

Section of sodium passage in C1 Section of sodium passage in C2 Sodium specific heat Reactivity induced by the control rods inserted into the corea Delayed neutron fraction Delayed neutron fraction Delayed neutron fraction Delayed neutron fraction Delayed neutron fraction Delayed neutron fraction Mean value of neutron production lifetimea Precursor decay constanta Precursor decay constanta Precursor decay constanta Precursor decay constanta Precursor decay constanta Precursor decay constanta

t2(T) A1 A2 CR(T)

gSCM b1 b2 b3 b4 b5 b6 G

g1 g2 g3 g4 g5 g6

a These parameters are usually denoted by lambda Greek letters instead of gamma Greek letters. In the present paper, gamma Greek letters are, however, used in order to avoid confusion with transition rates.

804

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Table 5 Initial conditions at time t ¼t0 for steady states of process variables. Process variable

Expression

Initial value

Unit

o(t0)

¼CM/K – ¼Te +R  w1  P(t0) ¼Te +R  w2  P(t0) ¼Te +(w1  P(t0))/ (2  A1  CR(T1(t0))  G(t0)) ¼Te +(w2  P(t0))/ (2  A2  CR(T2(t0))  G(t0)) – ¼(b1  P(t0))/(G  g1) ¼(b2  P(t0))/(G  g2) ¼(b3  P(t0))/(G  g3) ¼(b4  P(t0))/(G  g4) ¼(b5  P(t0))/(G  g5) ¼(b6  P(t0))/(G  g6)

6000 4000 1437.05 1773.15 873

rad s  1 kg m  2 s  1 K K K

G(t0) Tc,1(t0) Tc,2(t0) T1(t0) T2(t0) P(t0) C1(t0) C2(t0) C3(t0) C4(t0) C5(t0) C6(t0)

873

K

2,000,000,000.00 319,816,134.94 1,203,444,877.12 248,068,701.72 201,566,746.12 24,519,422.51 2,964,980.67

W W W W W W W

The sodium temperature in channel i, denoted Ti(t) [K], with i¼1, 2, is given by   d wi  PðtÞ 1 d 1 ðT ðtÞTe Þ  Ti ðtÞ ¼ GðtÞ þ ti ðTc,i ðtÞÞ i dt 2  Ai  ti ðTc,i ðtÞÞ  CR ðTi ðtÞÞ  GðtÞ GðtÞ dt

ð19Þ

m of data variables at time t is therefore y(t)¼ (Gm(t), Tm 1 (t), T2 (t), yCG(t), yCT(t), ySDL(t), Gd(t), Td1(t), Td2(t), a1(t), a2(t), Tc1(t), Tc2(t))T with yðtÞ A RL and L¼13. The data variables are modelled separately and defined by the following equations. The measurement results of the flow rate transmitter, denoted Gm(t) [kg m  2 s  1], are given by

Gm ðtÞ ¼ GðtÞ  I1 ðSSGðtÞ ¼ 1Þ þ  Gm  ðtÞ  I1 ðSSGðtÞ ¼ 0Þ þ Gmin  I1 ðSSGðtÞ ¼ 2Þ

ð22Þ with –Gm(t) which is the previous value of Gm(t), that is, – m G (t)¼Gm(t–e) with e which tends to 0 + (cf. Section 2.1). (For the numerical analyses using the formalism presented in Section 3.1, e ¼ Dt and –Gm(t+ Dt)¼Gm(t).) According to the components states given in Table 2, when the flow rate transmitter is in the mode of ‘‘normal’’ operation (i.e. SSG(t)¼1), it provides a correct result, that is, a measurement result equal to the sodium flow rate (i.e. Gm(t)¼G(t)); when it is in the ‘‘dangerous’’ failure mode (i.e. SSG(t)¼0), it always provides the previous measurement result (i.e. Gm(t)¼  Gm(t)); and when it is in the ‘‘safe’’ failure mode (i.e. SSG(t)¼2), it provides the threshold value (i.e. Gm(t)¼Gmin) or lower (i.e. Gm(t)oGmin which implies the same consequences on system behaviour as Gm(t)¼Gmin, according to the following equations, and as such is included in Gm(t)¼Gmin in Eq. (22)). The measurement results of the temperature transmitter i, denoted Tm i (t) [K], with i ¼1, 2, are given by

The power generated by the core, denoted P(t) [W], is given by ! 6 6 X d PðtÞ X PðtÞ ¼ gSCM  t  I1 ðSSCMðtÞ ¼ 5Þ bi þ gi  Ci ðtÞ dt G i¼1 i¼1

Tim ðtÞ ¼ ðTi ðtÞ þ Tic ðtÞÞ  I1 ðSSTiðtÞ ¼ 1,3, or 4ÞdDi ðtÞ  I1 ðSSTiðtÞ ¼ 3Þ þ dDi ðtÞ  I1 ðSSTiðtÞ ¼ 4Þ þ  Tim ðtÞ  I1 ðSSTiðtÞ ¼ 0Þ þTmax  I1 ðSSTiðtÞ ¼ 2Þ ð23Þ

ð20Þ

with the parameters of drift correction (data variables) Tci (t) defined hereafter, the drift deviation variables dDi(t) defined in m Section 4.5, and –Tm i (t) which is the previous value of Ti (t), that is, – m + Ti (t)¼Tm i (t–e) with e which tends to 0 (cf. Section 2.1). (For the numerical analyses using the formalism presented in Section 3.1, m e ¼ Dt and –Tm i (t + Dt) ¼Ti (t).) According to the components states given in Table 2 and the ‘‘intelligent transmitter’’ features, when the temperature transmitter i is in the mode of ‘‘normal’’ operation (i.e. SSTi(t)¼1), it provides a measurement result equal to the sodium temperature in channel i, plus the parameter of drift c correction (i.e. Tm i (t) ¼Ti(t)+ Ti (t)); when it is in a ‘‘degraded’’ mode of operation (i.e. SSTi(t) ¼3 or 4), the results are, in addition, c subject to negative or positive drifts (i.e. Tm i (t)¼Ti(t)+ Ti (t)–dDi(t) c or Tm i (t)¼ Ti(t)+ Ti (t) + dDi(t), respectively); when it is in the ‘‘dangerous’’ failure mode (i.e. SSTi(t) ¼0), it always provides the – m previous measurement result (i.e. Tm i (t)¼ Ti (t)); and when it is in the ‘‘safe’’ failure mode (i.e. SSTi(t)¼ 2), it provides the threshm old value (i.e. Tm i (t)¼Tmax) or greater (i.e. Ti (t)4Tmax which implies the same consequences on system behaviour as Tm i (t)¼ Tmax, according to the following equations, and as such is included in Tm i (t) ¼Tmax in Eq. (23)). The signals generated by the flow rate controller, denoted yCG(t), by the temperature controller, denoted yCT(t), and by the central controller, denoted ySDL(t), are given by

with the SCRAM state variable SSCM(t) defined in Section 4.2 (cf. Table 2). The precursor concentration i, denoted Ci(t) [W], with i¼1, 2, 3, 4, 5, 6, is given by d b C ðtÞ ¼ gi  Ci ðtÞ þ i PðtÞ dt i G

ð21Þ

Note that the introduction of an indicator function I1(A) has allowed the use of only one set of equations, which depend on components state variables as parameters. These equations also depend on deviation variables, and explicitly on time t (for Eq. (20)). Several derivatives of process variables are expressed according to other process variables and derivatives. In fact, it is sometimes not possible to provide explicit equations. However, these dependencies can be easily treated using the formalism presented in Section 3.1, when the variables are modelled separately, and by the help of the ‘‘meta-transitions’’ which make available the derivatives as additional variables. Finally Eqs. (16)–(21) are appropriate for the scenarios and analyses developed in the present paper, but are not applicable to any other cases. In particular, these equations apply only if the sodium temperatures in channels (T1(t) and T2(t)) are lower than 1156 K, and the fuel temperatures in channels (Tc,1(t) and Tc,2(t)) are lower than 3070 K. 4.4. Data variables Thirteen data variables are used: the measurement results of m the three transmitters (i.e. Gm(t), Tm 1 (t), and T2 (t)), the signals of the three controllers (i.e. yCG(t), yCT(t), and ySDL(t)), three stored data which correspond to measurement results at previous instants of time (i.e. Gd(t), Td1(t), and Td2(t)), two process data computed by the temperature transmitters to detect potential drifts (i.e. a1(t) and a2(t)), and two parameters of correction used to compensate for assumed drifts (i.e. Tc1(t) and Tc2(t)). The vector

yCG ðtÞ ¼ I1 ðGm ðtÞ r Gmin Þ  I1 ðSCGðtÞ ¼ 1Þ þI1 ðSCGðtÞ ¼ 2Þ

ð24Þ

yCT ðtÞ ¼ I1 ððT1m ðtÞ ZTmax Þ þI1 ðSCTðtÞ ¼ 2Þ

ð25Þ

or ðT2m ðtÞ ZTmax ÞÞ  I1 ðSCTðtÞ ¼ 1Þ

ySDL ðtÞ ¼ I1 ððyCG ðtÞ ¼ 1Þ or ðyCT ðtÞ ¼ 1ÞÞ  I1 ðSSDLðtÞ ¼ 1Þ þI1 ðSSDLðtÞ ¼ 2Þ

ð26Þ

According to Table 2, when the flow rate (resp. temperature) controller is in the mode of ‘‘normal’’ operation (i.e. SCG(t)¼1, resp. SCT(t) ¼1), it provides a correct signal, that is, equal to 1 if

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

the flow rate (resp. temperature) threshold is reached, and equal to 0 otherwise (i.e. yCG(t)¼I1(Gm(t)rGmin), resp. yCT(t)¼I1(Tm 1 (t) Z Tmax or Tm 2 (t)ZTmax)); when the central controller is in the mode of ‘‘normal’’ operation (i.e. SSDL(t)¼1), it provides a correct signal, that is, equal to 1 if at least one of the primary controllers sends a signal equal to 1, and equal to 0 otherwise (i.e. ySDL(t)¼ I1(yCG(t)¼1 or yCT(t) ¼1)); when a controller is in the ‘‘dangerous’’ failure mode (i.e. SCG(t)¼0, resp. SCT(t)¼0 and SSDL(t)¼ 0), it always provides a signal equal to 0 (i.e. yCG(t)¼0, resp. yCT(t) ¼0 and ySDL(t) ¼0), (according to the indicator function I1(A), these cases are implicit in Eqs. (24)–(26)) and therefore do not appear explicitly; and when a controller is in the ‘‘safe’’ failure mode (i.e. SCG(t) ¼2, resp. SCT(t) ¼2 and SSDL(t)¼2), it always provides a signal equal to 1 (i.e. yCG(t)¼1, resp. yCT(t)¼1 and ySDL(t)¼1). Eq. (23) shows that the results of the temperature transmitters may be subject to drifts (dD1(t) and dD2(t)) and, to compensate for these drifts, parameters of correction (Tc1(t) and Tc2(t)) are added. These parameters of drift correction are computed using the following ‘‘intelligent transmitter’’ features: i. The three transmitters (SG, ST1, and ST2) store their previous measurement results, denoted Gd(t), Td1(t), and Td2(t), respectively. ii. The flow rate transmitter (SG) transmits its current and previous measurement results (Gm(t) and Gd(t)) to the temperature transmitters. iii. Each temperature transmitter i (STi), with i¼ 1, 2, computes process data, denoted ai(t), according to its current and d previous measurement results (Tm i (t) and Ti (t)), and the data m received from the flow rate transmitter (G (t) and Gd(t)), used to detect potential drifts in its measurement results. iv. The temperature transmitters (ST1 and ST2) transmit to each other their process data (a1(t) and a2(t)), and their current m measurement results (Tm 1 (t) and T2 (t)). v. Each temperature transmitter i (STi), with i¼ 1, 2, computes a parameter of correction, denoted Tci (t), according to the process data (a1(t) and a2(t)), and the current measurement m results (Tm 1 (t) and T2 (t)), used to compensate for assumed drifts in its measurement results. The stored data that correspond to the previous measurement results of the flow rate transmitter, denoted Gd(t) [kg m  2 s  1], and of the temperature transmitters i, denoted Tdi (t) [K], with i¼1, 2, are given by Gd ðtÞ ¼  Gm ðtÞ

ð27Þ

Tid ðtÞ ¼  Tim ðtÞ

ð28Þ

(For the numerical analyses using the formalism presented in Section 3.1, Gd(t + Dt) ¼Gm(t) and Tdi (t + Dt) ¼Tm i (t).) Note that it is also possible to define Gd(t) and Tdi (t) using a larger number of

805

prior values of Gm(t) and Tm i (t), respectively, using conditions on time t (cf. for example the function given in Appendix A), which can be easily treated using the formalism presented in Section 3.1. The process data ai(t), with i¼1, 2, evaluate the derivative of the measurement results of the temperature transmitter i, accordd ing to Tm i (t) and Ti (t), and compare it with the theoretical m d derivative obtained by Eq. (19), using Tm i (t), Ti (t), G (t), and Gd(t), instead of the ‘‘real’’ (and unknown) values of the process variables. To this end, the parameters given in Table 4 are assumed to be known, as well as the initial conditions given in Table 5. Moreover, according to Eqs. (18), (20), and (21), and Table 5, it is shown that the process variables Tc,i(t), with i¼ 1, 2, and P(t), are constant and equal to Tc,i(t0) and P(t0), respectively (and independently of the deviation variables), up to the SCRAM activation (the period of time after the SCRAM activation is not relevant for the following analyses), (cf. proof in Appendix B). All the parameters are therefore available to compute ai(t), with i¼1, 2:

ai ðtÞ ¼

Tim ðtÞTid ðtÞ

e

Bi ðtÞ

ð29Þ

with Bi ðtÞ ¼

wi  Pðt0 Þ 2  Ai  ti ðTc,i ðt0 ÞÞ  CR ðTim ðtÞÞ  Gm ðtÞ  m  G ðtÞGd ðtÞ 1 ðT m ðtÞTe Þ  þ Gm ðtÞ  e ti ðTc,i ðt0 ÞÞ i

ð30Þ

with e which tends to 0 + . (For the numerical analyses using the formalism presented in Section 3.1, and according to Eqs. (27) and (28), e ¼ Dt.) Analyses have then shown that measurement results m m Tm 1 (t), T2 (t), and G (t), can be further characterised using process data a1(t) and a2(t), as shown in Table 6. According to Eq. (19) and Table 5, the temperatures T1(t) and T2(t) are almost equal (and independently of the deviation variables), up to the SCRAM activation (cf. empirical results obtained for the scenarios given in Section 5.3). It is then

Table 6 Relationships between process data ai(t), and measurement results Tm i (t) and Gm(t), with i¼ 1, 2. m Measurement resultsa Tm i (t) and G (t)

Process data ai(t)

m Tm i (t) ¼ Ti(t) and G (t)¼ G(t) m Tm i (t) o Ti(t) or G (t) 4G(t) m Tm i (t) 4Ti(t) or G (t) o G(t)

ai(t) ¼0 ai(t) o0 ai(t) 40

a Other second order combinations of events such as ‘‘Tm i (t) oTi(t) and m Gm(t)o G(t)’’ or ‘‘Tm i (t) 4Ti(t) and G (t) 4G(t)’’ yield indeterminate results for ai(t) because they depend on the degrees of difference between Tm i (t) and Ti(t), and between Gm(t) and G(t).

Table 7 Definition of parameters of drift correction Tci (t) [K], according to ai(t), with i¼1, 2. Criteriaa on ai(t)

Accepted hypothesisb

Correction parameters Tci (t)

a1(t) o–aref and a2(t) o –aref a1(t) 4 aref and a2(t) 4 aref 9a1(t)94 aref and 9a2(t)9o aref 9a1(t)9o aref and 9a2(t)94 aref a1(t) o– aref and a2(t) 4 aref a1(t) 4 aref and a2(t) o –aref

Gm(t) 4G(t) Gm(t) oG(t) m Tm 1 (t) o T1(t) or T1 (t) 4T1(t) m Tm 2 (t) o T2(t) or T2 (t) 4T2(t) m Tm 1 (t) o T1(t) and T2 (t) 4T2(t) m Tm 1 (t) 4T1(t) and T2 (t)o T2(t) m Tm 1 (t) ¼T1(t) and T2 (t) ¼T2(t)

Tci (t) ¼  Tci (t) for i¼ 1, 2 Tci (t) ¼  Tci (t) for i¼ 1, 2  c m c Tc1(t) ¼ Tm 2 (t)  T1 (t) and T2(t)¼ T2(t) m Tc1(t) ¼  Tc1(t) and Tc2(t) ¼Tm 1 (t)  T2 (t) m c m m Tc1(t) ¼ ((Tm 2 (t)  T1 (t))/2) and T2(t) ¼ ((T1 (t)  T2 (t))/2) m c m m Tc1(t) ¼ ((Tm 2 (t)–T1 (t))/2) and T2(t) ¼((T1 (t)–T2 (t))/2) Tci (t) ¼  Tci (t) for i¼ 1, 2

Otherwise a

With the parameter aref discussed in Section 5.4. The accepted hypotheses are defined according to Table 6, based on a priori equiprobability between the following events: Gm(t) o G(t), Gm(t) 4G(t), Tm i (t) oTi(t), and Tm i (t) 4Ti(t) for i ¼1, 2. b

806

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

proposed to use the process data a1(t) and a2(t) to detect potential m drifts in the measurement results Tm 1 (t) and T2 (t) and, if relevant, to define the parameters of drift correction Tc1(t) and Tc2(t) in such a way that both measurement results are equal. According to a1(t), a2(t), and Table 6, different hypotheses can be assumed with regard to drifts in measurement results, and they are used to define the parameters of drift correction Tc1(t) and Tc2(t), according to the rules given in Table 7. If drifts are assumed only in the measurement results of the temperature transmitter 1 (resp. transmitter 2), then the parameter of drift correction Tc1(t) (resp. m Tc2(t)) is changed in such a way that Tm 1 (t) becomes equal to T2 (t), c m m c m m that is, T1(t)¼T2 (t)–T1 (t) (resp. T2(t)¼ T1 (t)–T2 (t)). If drifts are assumed in the measurement results of both temperature transmitters 1 and 2, then the drift corrections Tc1(t) and Tc2(t) are both m changed in such a way that Tm 1 (t) and T2 (t) become equal to m c m m c m (Tm 1 (t)+ T2 (t))/2, that is, T1(t)¼ ((T2 (t)–T1 (t))/2) and T2(t)¼((T1 (t)– m T2 (t))/2). In the other cases, the drift corrections are unchanged (i.e. Tc1(t) ¼ –Tc1(t) and Tc2(t)¼ –Tc2(t)). The parallel between the data variables defined in this section and the model for digital-based transmitters shown in Fig. 4 may be drawn. For example, assuming the modelling of the temperature transmitter i (STi), the measurements results Tm i (t) correspond to the subset (M) of data variables, the stored data Tdi (t) to the subset (D), the process data ai(t) to the subset (P), the parameters of correction Tci (t) to the subset (C), and the transmitted information to the other system components are comprised of Tm i (t) and ai(t), which correspond to the subset (I). The parameter aref (used in the criteria to define the parameters of drift correction Tci (t), cf. Table 7) can be assumed to be either a fixed parameter, that is, not represented as a variable in the modelling framework, or an input variable which, for example, can be changed by a user, and therefore is assumed to be a data variable received from other system components (i.e. subset \(T) of vector y(t)). In further development, additional diagnostic information could be assumed, for example to detect a ‘‘non-operating’’ mode of the transmitter if the parameter of drift correction Tci (t) reaches a certain level. Similarly no parameter of reconfiguration (subset (R)) is used for the present case study. Such parameters could, for example, refer to the command of a cold stand-by if the transmitter is detected as ‘‘nonoperating.’’

4.5. Deviation variables Three deviation variables are used: the pump torque deviations (i.e. dM(t), cf. Eq. (16)), and the temperature transmitter drifts (i.e. dD1(t) and dD2(t), cf. Eq. (23)). The vector of deviation variables at time t is therefore e(t) ¼(dM(t), dD1(t), dD2(t))T, with eðtÞ A RQ and Q¼3. The deviation variables are modelled separately and defined by the following stochastic processes. The pump torque deviation, denoted dM(t) [N m], is defined as

dM ðt þ dtÞ ¼ dM ðtÞþ EM ðaM  dt, bM Þ  I1 ðdM ðtÞo CM Þ  I1 ðSPMðtÞ ¼ 3Þ ð31Þ with the pump state variable SPM(t) defined in Section 4.2 (cf. Table 2), and the random variable EM(aM  dt, bM) which follows a Gamma distribution of shape parameter aM  dt¼ 0.5  dt and scale parameter 1/bM ¼400.0 [N m]. The parameter CM is the nominal pump torque given in Table 4. According to Eq. (16), the total pump torque is therefore equal to CM–dM(t) when SPM(t)¼3. Since the increment of the pump torque deviation in time interval [t, t + dt], defined by EM(aM  dt, bM)  I1(dM(t)oCM)  I1(SPM(t)¼3), is independent on time t (but, of course, depends on dt), this stochastic process is stationary, given I1(dM(t)oCM)  I1(SPM(t)¼3). (For the

numerical analyses using the formalism presented in Section 3.1, dt¼ Dt.) The temperature transmitter drifts, denoted dDi(t) [K], with i¼1, 2, are defined as

dDi ðt þ dtÞ ¼ dDi ðtÞ þ ED ðaD  dt, bD ðdDi ðtÞ,Ti ðtÞÞÞ  I1 ðSSTiðtÞ ¼ 3 or 4Þ ð32Þ with the temperature transmitter state variable SSTi(t) defined in Section 4.2 (cf. Table 2), the temperature process variable Ti(t) defined in Section 4.3, and the random variable ED(aD  dt, bD(dDi(t), Ti(t))) which follows a Gamma distribution of shape parameter aD  dt¼ 1.0  dt and scale parameter 1/bD(dDi(t), Ti(t))¼ (Ti(t)  (dDi(t)+20.0))/3.0  104 [K], with i¼ 1, 2. (For the numerical analyses using the formalism presented in Section 3.1, dt¼ Dt.)

5. Case study: modelling and analyses 5.1. Dynamic reliability modelling The system has been modelled using the formalism presented in Section 3.1, and the corresponding Petri net is shown in Fig. 6. In accordance with this formalism, each variable has been modelled separately (except for the precursor concentrations process variables (Ci(t), with i¼1, y, 6) which have been modelled as one vector, denoted C(t)). The process variables are shown on the left-hand side of Fig. 6, and can be distinguished by the places in plain grey circles, that are associated with ‘‘metatransitions,’’ (cf. caption in Fig. 1). The ‘‘meta-transition’’ associated to the flow rate variable (G(t)) has been explicitly drawn because the derivative included in this ‘‘meta-transition’’ (G0 (t)) is used to compute the evolution of the sodium temperatures (T1(t) and T2(t)). The deviation variables can be distinguished by the plain white circles, the components state variables by the dashed white circles, and the data variables by the plain light grey circles that are associated to ‘‘normal’’ transitions (cf. Fig. 1). Finally the time t variable has been also depicted and directly affects the SCRAM state variable (SSCM(t), cf. Table 3) and the power process variable (P(t), cf. Eq. (20)). To improve the legibility of Fig. 6, the specifications of the transition guards have not been drawn. However, it is assumed that each transition is fired at every time step Dt, and has therefore a guard denoted si[Dt], with the time instants si which respect the specific order given in Section 3.1: the evolutions of process and deviation variables are computed first, then the components states variables, the time t variable, the process and deviation variables, and finally the data variables. In addition, among the variables of the same type, the derivative of the flow rate variable (G0 (t)) has to be computed before the derivatives of the sodium temperatures (T1(t) and T2(t)), because the latter depend on the first, and the data variables should, for example, be computed according to this order: Gd(t), Td1(t), Td2(t), Gm(t), Tm 1 (t), Tm 2 (t) (because the measurement results have to be stored according to Eqs. (27) and (28), before computing the new results according to Eqs. (22) and (23)), a1(t), a2(t), Tc1(t), Tc2(t), yCG(t), yCT(t), ySDL(t). Note that the parameters of drift corrections computed at time t (Tc1(t) and Tc2(t)) are therefore applied to the measurement results at time t + Dt. The specifications of the output arcs have not been drawn either. According to Sections 2.2 and 3.1, the output arc that corresponds to each process variable x(t)k (i.e. kth component of vector x(t)) should be specified by the right-hand side of the expression x(t + Dt)k ¼x(t)k + Dt  x0 (t)k, with the derivatives x0 (t)k given by the respective Eqs. (16)–(21). Similarly the output arc

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

807

Fig. 6. Petri net for the dynamic reliability modelling of the case study.

that corresponds to each deviation variable e(t)k should be specified by the right-hand side of the expression e(t+ Dt)k ¼ e(t)k + dE(t)k, with the random increment dE(t)k given by the respective Eqs. (31) and (32). The output arc that is used to change the value of each component state variable i(t)k should be specified by the righthand side of the expression i(t + Dt)k ¼P[ik],i(t)(x(t), y(t), e(t), t)T  [ik], with [ik] defined by all the possible values of variable i(t)k, as given in Table 2, and P[ik],i(t)(x(t), y(t), e(t), t) the random vector defined by the procedure described in Section 2.2, using the transition rates given in Table 3. For example, the flow rate transmitter state variable (SSG(t)), (for this example, k¼2 and i(t)k ¼SSG(t) according to Section 4.2), has ESSG ¼3 possible states, and [SSG] ¼(1, 0, 2)T, according to Table 2. The random vector P[SSG],SSG(t), (for this example, the variables x(t), y(t), e(t), and t, are not relevant), is defined according to the value of SSG(t), Table 3, and the procedure described in Section 2.2: if SSG(t)¼1 then q[SSG],SSG(t) ¼((1–Dt  (lSG,0 + lSG,2)), Dt  lSG,0, Dt  lSG,2)T, if SSG(t) ¼0 then q[SSG],SSG(t) ¼(0, 1, 0)T, and if SSG(t)¼2 then q[SSG],SSG(t) ¼(0, 0, 1)T, then R[SSG],SSG(t) is a discrete random variable with sample space {1, 2, 3} and probability mass function defined by q[SSG],SSG(t), finally P[SSG],SSG(t) is a vector with all its components equal to 0, except for its component indexed by the realisation of R[SSG],SSG(t), that is equal to 1. The output arc that corresponds to each data variable y(t)k should be directly specified by the respective Eqs. (22)–(30), and Table 7 for the parameters of drift correction (Tc1(t) and Tc2(t)), with  m m e ¼ Dt and  Gm(t+ Dt)¼Gm(t),  Tm T2 (t+ Dt)¼Tm 1 (t+ Dt)¼T1 (t), 2 (t), according to the formalism presented in Section 3.1. Finally for the output arc that corresponds to the time t variable, the expression of the output arc is simply t + Dt.

Following the model shown in Fig. 6, the system has been modelled and its evolution simulated, using CPN Tools [86,87], a software developed by the University of Aarhus for editing, simulating and analysing timed and untimed coloured Petri nets. This tool possesses interesting properties which governed our selection:

 Ability to represent coloured properties.  Expressive semantics for variable and expression definitions 

   

(notably for the specification of the output arcs), including stochastic properties. Ability to model hierarchical Petri Nets by assigning a separate Petri Net to a transition (useful to represent ‘‘meta-transitions’’) or by creating fusion places (which may improve the clarity of the model). Supports analyses by Monte Carlo simulations. Editing capabilities and ability to define monitoring functions (for example, to monitor the evolution of any system variable). User friendly interface, including a graphical editor with embedded syntax chequer, and contextual error messages. Free computer tool. Some drawbacks can, however, be identified:

 The tool does not handle real numbers and those have therefore to be ‘‘encoded’’ using integers.

 The simulations are often relatively time-consuming. Several analyses have been performed and are described in the following sections. For all these analyses, the initial time is t0 ¼0.0 s, and the time step is set to 1.0 s, that is Dt ¼1.0 s.

808

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Fig. 7. Examples of pump torque (CM  dM(t)) [N m] evolution, when the pump is in a ‘‘degraded’’ mode of operation.

Fig. 9. Evolutions of the sodium flow rate (G(t)) [kg m  2 s  1] in the scenarios ‘‘a’’ and ‘‘b’’.

Fig. 10. Evolutions of the sodium temperatures (T1(t) and T2(t)) [K] in the scenarios ‘‘a’’ and ‘‘b’’. Fig. 8. Examples of temperature transmitter drifts (dDi(t)) [K] evolution, when the temperature transmitter is in a ‘‘degraded’’ mode of operation.

5.2. Examples of deviation variables evolution Examples of stochastic evolution of the pump torque (CM–

dM(t)) are shown in Fig. 7, according to the pump torque deviation (dM(t)) expressed by Eq. (31), starting at time t ¼0 s with the pump in the ‘‘degraded’’ mode of operation and staying in this state (SPM(t)¼ 3 for t Z0). Examples of stochastic evolution of the temperature transmitter drifts (dDi(t)) are shown in Fig. 8, according to Eq. (32), starting at time t ¼0 s with the corresponding temperature transmitter that is subject to drifts and staying in this state (STTi(t)¼3 or 4 for t Z0), and with the corresponding sodium temperature (Ti(t)) that is assumed to increase linearly from Ti(0 s) ¼875 K to Ti(60 s) ¼935 K.

5.3. Examples of scenarios At time t ¼t0 ¼0.0 s, the system components are in the ‘‘normal’’ modes of operation (all the components state variables are equal to 1) except for those specified hereafter; the process variables are in the initial conditions defined by Table 5; the data variables for the measurement results and signals are defined by Eqs. (22)–(26), and the other data variables are equal to 0.0; and the deviations variables are equal to 0.0. The following scenarios have then been simulated by setting the evolution of components state variables, and the corresponding evolutions of the sodium

flow rate (G(t)) and temperatures (T1(t) and T2(t)) are shown in Figs. 9–12: a. The pump stays in the ‘‘normal’’ mode of operation up to time t¼30 s (SPM(t) ¼1 for t0 rt o30 s) where a full failure of the pump occurs (SPM(t)¼0 for tZ30 s), and the SCRAM is activated at time t ¼36 s (SSCM(t)¼1 for t0 rt o36 s and SSCM(t)¼5 for t Z36 s) with the corresponding process variables (obtained by simulations): T1(36 s) ¼924.5 K, T2(36 s)¼ 922.5 K, and G(36 s) ¼3305.9 kg m  2 s  1. b. The pump stays in the ‘‘normal’’ mode of operation up to time t¼30 s (SPM(t) ¼1 for t0 rt o30 s) where a full failure of the pump occurs (SPM(t)¼0 for t Z30 s), and the SCRAM is never activated (SSCM(t) a5 for t0 rt). c. The pump is in the ‘‘degraded’’ mode of operation at time t0 and stays in this state (SPM(t) ¼3 for t Zt0), and the SCRAM is activated at time t ¼45 s (SSCM(t)¼1 for t0 rt o45 s and SSCM(t)¼5 for t Z45 s) with the corresponding process variables (obtained by simulations): T1(45 s) ¼923.5 K, T2(45 s) ¼ 923.1 K, and G(45 s) ¼3404.4 kg m  2 s  1. d. The pump is in the ‘‘degraded’’ mode of operation at time t0 and stays in this state (SPM(t) ¼3 for t Zt0), and the SCRAM is activated at time t ¼68 s (SSCM(t)¼1 for t0 rt o68 s and SSCM(t)¼5 for t Z68 s) with the corresponding process variables (obtained by simulations): T1(68 s) ¼930.1 K, T2(68 s) ¼ 929.6 K, and G(68 s) ¼3342.6 kg m  2 s  1. e. The pump is in the ‘‘degraded’’ mode of operation at time t0 and stays in this state (SPM(t) ¼3 for t Zt0), and the SCRAM is never activated (SSCM(t)a5 for t0 rt). The evolutions of the sodium flow rate (G(t)) in the scenarios ‘‘a’’ and ‘‘b’’ are, respectively, denoted G(t) . a and G(t) . b, and are

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

809

difference between T1(t) and T2(t) globally increases as a function of time, but this difference is less than 1.00 K for, at least, t r90 s). 5.4. Reliability analyses

Fig. 11. Evolutions of the sodium flow rate (G(t)) [kg m  2 s  1] in the scenarios ‘‘c,’’ ‘‘d,’’ and ‘‘e’’.

Fig. 12. Evolutions of the sodium temperatures (T1(t) and T2(t)) [K] in the scenarios ‘‘c,’’ ‘‘d,’’ and ‘‘e’’.

shown in Fig. 9. According to the steady state conditions, the flow rate is constant up to the full failure of the pump at time t ¼30 s, and then decreases according to Eqs. (16) and (17). Note that, according to these equations, the SCRAM activation has no effect on the sodium flow rate, thus G(t) . a¼ G(t) . b. The evolutions of the sodium temperatures (T1(t) and T2(t)) in the scenarios ‘‘a’’ and ‘‘b’’ are, respectively, denoted T1(t) . a, T2(t) . a, T1(t) . b, and T2(t) . b, and are shown in Fig. 10. According to the steady state conditions, the temperatures are constant up to the full failure of the pump at time t¼30 s, then increase according to Eqs. (16)–(21) and, once the SCRAM is activated at time t¼ 36 s (for the scenario "a"), decrease after a certain delay. Note that the sodium temperatures are almost equal before the SCRAM activation (the difference between T1(t) and T2(t) is less than 0.02 K before the full failure of the pump, and then this difference increases). Once the SCRAM is activated, T1(t) decreases faster than T2(t). The evolutions of the sodium flow rate (G(t)) in the scenarios ‘‘c,’’ ‘‘d,’’ and ‘‘e’’ are, respectively, denoted G(t) . c, G(t) . d, and G(t) . e, and are shown in Fig. 11. According to Eqs. (16) and 17, and the stochastic pump torque deviation defined by Eq. (31), the flow rate decreases randomly in these three scenarios and therefore they do not yield the same values. The evolutions of the sodium temperatures (T1(t) and T2(t)) in the scenarios ‘‘c,’’ ‘‘d,’’ and ‘‘e’’ are respectively denoted T1(t) . c, T2(t) . c, T1(t) . d, T2(t) . d, T1(t) . e, and T2(t) . e, and are shown in Fig. 12. According to Eqs. (16)–(21), and the stochastic pump torque deviation defined by Eq. (31), the temperatures increase randomly in these three scenarios. Note that the sodium temperatures are almost equal before the SCRAM activation (the

To perform reliability analyses, a safety domain D has to be defined. For this application, it logically depends on the sodium temperatures in the channels (T1(t) and T2(t)), and the safety boundary is defined by TD ¼933 K. However, according to the scenarios ‘‘a’’ and ‘‘b’’ presented in Section 5.3, it is shown that, in case of a full failure of the pump, the sodium temperatures increase so quickly that, even if the safety function is well performed, the boundary TD may be reached for a couple of seconds. The definition of the safety domain D is therefore completed by a time condition and a ‘‘hazardous event’’ is then assumed to occur if and only if at least one of the sodium temperatures (T1(t) or T2(t)) is equal to or greater than TD ¼933 K for more than 5 s. The safety function has to be defined according to the threshold values of sodium temperature (Tmax), and flow rate (Gmin), which should activate the SCRAM. For this application, the values Tmax ¼ 923 K and Gmin ¼3345 kg m  2 s  1 are used. For example, in the scenarios ‘‘a’’ and ‘‘b’’ presented in Section 5.3, once a full failure of the pump has occurred, both threshold values are reached almost at the same time, and after about 6 s. In the scenarios ‘‘c,’’ ‘‘d,’’ and ‘‘e’’ it is shown that, when the pump is in the ‘‘degraded’’ mode of operation, the temperature threshold is reached first (at a random time, according to the pump torque deviation), and then the flow rate threshold is reached approximately 2–5 s later. In these examples, the scenario ‘‘c’’ corresponds to a case where the SCRAM is activated as soon as the temperature threshold is reached, and the scenario ‘‘d’’ to a case where the SCRAM is activated only when the flow rate threshold is reached (for example, in case of both temperature transmitters are in ‘‘dangerous’’ failure modes). Since the system components have ‘‘safe’’ failure modes (and the temperature transmitters may be subject to ‘‘safe’’ drifts), it is also possible that the threshold values are ‘‘detected’’ by components, even if the temperatures and the flow rate are, actually, still under the threshold values. The SCRAM may therefore be incorrectly activated, which implies a ‘‘spurious trip.’’ Since the temperature transmitters may be subject to drifts, it is however probably inappropriate to count as a ‘‘spurious trip’’ a scenario where the SCRAM is activated just a few degrees lower than the temperature threshold. For this application, a ‘‘spurious trip’’ is therefore assumed if and only if the SCRAM is activated while both sodium temperatures are lower than Tmax  5 K¼918 K. When neither a ‘‘hazardous event’’ nor a ‘‘spurious trip’’ has occurred, then the system is in its useful life. During this period, if the SCRAM is activated while at least one sodium temperature is greater than Tmax–5 K¼918 K, but none of the sodium temperatures have exceeded TD ¼933 K for more than 5 s, then an ‘‘activation under control’’ is assumed. To sum up, the three possible scenarios (‘‘hazardous event,’’ ‘‘spurious trip,’’ and ‘‘activation under control’’) are classified according to the criteria shown in Table 8. The reliability analyses are performed according to an initiating event that is the occurrence of a partial failure (the ‘‘degraded’’ mode of operation) of the pump (SPM(t)¼3). That is, at time t¼ t0, the system components are in the ‘‘normal’’ modes of operation (all the components state variables are equal to 1) except for the pump which is in the ‘‘degraded’’ mode of operation (SPM(t0) ¼3); the process variables are in the initial conditions defined by Table 5; the data variables for the measurement results and signals are defined by Eqs. (22)–(26), and the other data variables are equal to 0.0; and the deviations variables are equal to 0.0.

810

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Table 8 Classification of the scenarios. Scenario

Criteria

Hazardous event Spurious trip

T1(t) ZTD or T2(t) ZTD for more than 5 s with TD ¼ 933 K

Activation under control

T1(t) oTmax–5 K and T2(t) oTmax  5 K while SCRAM becomes activated, with Tmax ¼ 923 K Other cases i.e. T1(t) Z Tmax–5 K or T2(t) ZTmax–5 K while SCRAM becomes activated, and neither T1(t) Z TD nor T2(t) ZTD for more than 5 s

is extended, as shown by Table 10, and ‘‘hazardous events’’ are automatically more likely due to the fact that probabilities of ‘‘dangerous’’ failure modes increase with time. (In addition, note that the SCRAM full failure, which always implies a ‘‘hazardous event,’’ occurs according to a transition rate that also increases with time.) These observations are shown by Figs. 13 and 14 which show the probabilities that a ‘‘hazardous event’’ and a ‘‘spurious strip,’’ respectively, has occurred before time t, and Fig. 15 which shows the probability that the system is still in its useful life at time t, without and with ‘‘intelligent transmitter’’ features (with aref ¼ + N and aref ¼1.0, respectively).

Table 9 Results of reliability analyses: percentagea of each scenario. Scenario Hazardous event (7 0.09)b Spurious trip (7 0.26)b Activation under control (7 0.26)b

aref ¼ + N

aref ¼ 5.0

aref ¼3.0

aref ¼ 1.0

aref ¼ 0.5

2.4

3.4

3.2

3.4

3.2

51.5

47.8

44.5

44.0

44.2

46.1

48.8

52.3

52.7

52.5

a

Obtained through 100,000 Monte Carlo simulation trials for each column. 90% confidence intervals, around the averages given in this table, assuming Binomial distributions approximated by Normal distributions. b

Fig. 13. Probability that a ‘‘hazardous event’’ has occurred before time t, without and with ‘‘intelligent transmitter’’ features (with aref ¼ + N and aref ¼ 1.0, respectively).

Table 10 Results of reliability analyses: average useful life.a

Average useful life (70.07)b

aref ¼ + N

aref ¼ 5.0

aref ¼3.0

aref ¼ 1.0

aref ¼ 0.5

48.86 s

51.32 s

51.57 s

51.63 s

51.57 s

a

Obtained through 100,000 Monte Carlo simulation trials for each column. 90% confidence intervals, around the averages given in this table, assuming Normal distributions. b

To assess the effects of the ‘‘intelligent transmitter’’ features, several cases are analysed according to the value of aref (cf. Section 4.4 and Table 7). When aref ¼ +N, the parameters of drift correction Tc1(t) and Tc2(t) are always equal to 0.0, according to their values at time t ¼t0 and Table 7, that is, no drift correction is assumed. Otherwise, the smaller aref and the more the parameters of drift correction are subject to changes, according to the rules defined in Table 7. The analyses have been performed using Monte Carlo simulations, for aref ¼ + N, aref ¼5.0, aref ¼3.0, aref ¼1.0, and aref ¼0.5 (100,000 trials were used for each of these cases). The respective percentages obtained for each scenario (‘‘hazardous event,’’ ‘‘spurious trip,’’ and ‘‘activation under control’’) are shown in Table 9, and the respective average useful lives are shown in Table 10. Table 9 shows that the use of the ‘‘intelligent transmitter’’ features allow a greater percentage of ‘‘activation under control,’’ by reducing the percentage of ‘‘spurious trip.’’ In fact, the drift corrections compensate for ‘‘safe’’ (positive) drifts of the temperature transmitters which are significant ‘‘spurious trip’’ contributors. On the other hand, the percentage of ‘‘hazardous event’’ is not reduced by the drift corrections, notably because even in case of ‘‘dangerous’’ (negative) drifts of the temperature transmitters, the SCRAM can still be activated by the detection of a low flow rate (the flow rate monitoring acts as a redundancy for the safety function). It can even be shown that the use of the drift corrections slightly increases the percentage of ‘‘hazardous event.’’ This is a direct consequence of the reduction of ‘‘spurious trips.’’ By preventing ‘‘spurious trips,’’ the useful life of the system

Fig. 14. Probabilities that a ‘‘spurious strip’’ has occurred before time t, without and with ‘‘intelligent transmitter’’ features (with aref ¼ + N and aref ¼ 1.0, respectively).

Fig. 15. Probability that the system is still in its useful life at time t, without and with ‘‘intelligent transmitter’’ features (with aref ¼ + N and aref ¼ 1.0, respectively).

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

Using the drift corrections, the choice of the parameter aref then affects the system performances according to non-monotonic relationships. With regard to the percentages of ‘‘activation under control’’ and ‘‘spurious trip,’’ and to the average useful life, an optimal value of aref can be found around 1.0, according to the cases shown in Table 9. In fact, when the parameter aref is smaller, the drift corrections are more sensitive and the probability of accepting a wrong hypothesis from the rules defined in Table 7 increases. In particular, accurate measurement results can then be biased by parameters of drift corrections that are defined according to inaccurate measurement results (notably in case of multiple failures, for example when both a temperature transmitter and the flow rate transmitter are not in ‘‘normal’’ modes of operation).

6. Conclusion To deal with the specific characteristics of digital-based transmitters which are able to process, store, exchange data, and perform advanced functionalities such as self-diagnosis and error measurement correction, the mathematical framework of the dynamic reliability has been extended by introducing data variables. In addition, deviation variables have been used to model system components errors such as drifts. A formalized Petri net approach has then been proposed in order to model, simulate, and analyse the dynamic reliability of such systems. To this end, numerical methods are used, based on Taylor series and finite differences. Following this formalism, it has been notably possible to propose a versatile model for digital-based transmitters, and to model and analyse a substantial case study. The case study is comprised of 8 components state variables, for a total number of 36  52 ¼18,225 possible states, according to Table 2; 13 process variables, although only two of them are directly affected by the other types of variables (the angular speed of the pump and the power generated by the core) and could therefore theoretically be used to exhaustively express all the process dynamics; 13 data variables, including three binary and ten continuous variables; and 3 deviation variables. Assuming a discretization of each continuous variable by a number of m states, the total number of states that describes the complete system is 36  52  m2  23  m10  m3, at each time t. Even for a small value of m, for example m¼16 (which is certainly not enough to provide a correct analysis of the system), the total number of states of the complete system reaches 1.68  1023, that is, more than the number of grains of sand in all the world’s beaches and deserts, and more than the number of stars in the known universe [88]. Even faced with these very large numbers of states, the complete system has been modelled using Petri nets with, according to the proposed formalism, a number of places and transitions that are both equal to the number of system variables (including time t), that is 38, plus the number of places and transitions used in the ‘‘meta-transitions’’ for process and deviation variables, that is 16, for a total of 54. Because the modelling size is linearly dependent on the number of variables, there is no real scale limitation to the modelling of more complex systems. However, the time required to perform the simulations and the analyses is another concern. In fact, even though a software like CPN Tools is flexible enough to meet the main formalism requirements, it is not optimised for these simulation tasks, which are therefore relatively time-consuming. For the developed case study, the processing time required to simulate one scenario using this software is about 4 s with a 2.20 GHz processor. However, empirical tests with experimental tools have shown that this processing time could be reduced by factors of 102–103 using optimised software dedicated

811

to these analyses. The development of such a software tool would then be an important task to extend the use of the proposed approach. Finally it is therefore concluded that the presented framework provides a way of modelling and analysing complex systems, including digital-based transmitters, able to deal efficiently with large-scale problems, while providing modelling flexibilities and a visual interface. Such criteria for an integrated dynamic reliability platform in PRA have been discussed by Labeau et al. [35]. It is then believed that the proposed formalism likely provides good prospects for such a tool. Further developments would, however, be of interest as the improvement of the Monte Carlo methods, for example using biasing techniques and adaptive time steps. By reducing the simulation time it would then be practical to perform uncertainty analyses. Other prospects for developments include the integration of related works, notably on software reliability, digital communication dependability, and human operations. Finally formal procedures should be established in order to integrate these frameworks in PRA.

Acknowledgments This work was performed at The Ohio State University (OSU) as part of a collaboration between the Department of Mechanical Engineering (Nuclear Engineering Program) of OSU, the French National Institute for Industrial Environment and Risk (INERIS), and the Troyes University of Technology (UTT). During this work, Florent Brissaud was supported by the French Ministry of Higher Education and Research, and the French Ministry of Ecology, Energy, and Sustainable Development.

Appendix A. Example of function used to store data variables The following function, denoted Store(y(r), –y(D), t), may be used to store K values of y(t)(r), with K A N, which correspond to the time instants tr (k–1)  y, with y A R þ and k¼1, 2, y, K, in the vector y(t)(D), with y(t)(D),1 ¼tr and y(t)(D),k + 1 ¼y(tr–(k–1)  y)(r), (y(t)(D),k + 1 is (k+1)th component of vector y(t)(D)): yðtÞðDÞ ¼ ðtr , yðtr ÞðrÞ , yðtr yÞðrÞ ,. . ., yðtr ðK1Þ  yÞðrÞ Þ

¼ StoreðyðrÞ , 2 yðDÞ , tÞ 2 ¼ ðt, yðtÞðrÞ , yðtÞðDÞ,2 , . . .,

2

yðtÞðDÞ,K Þ

if tr þ y r t

2

and yðtÞðDÞ otherwise

ðA:1Þ

Appendix B. Evolution of process variables According to Table 5, the fuel temperature in channel i at time t¼ t0, denoted Tc,i(t0) [K] with i¼1, 2, is Tc,i(t0)¼Te +R  wi  P(t0), and the precursor concentration i at time t ¼t0, denoted Ci(t0) [W] with i¼1, 2, 3, 4, 5, 6, is Ci(t0) ¼(bi  P(t0))/(G  gi). According to Eq. (21), the derivative of the precursor concentration i at time t ¼t0, with i¼1, 2, 3, 4, 5, 6, is therefore, d b bi b C ðt0 Þ ¼ gi  Ci ðt0 Þ þ i Pðt0 Þ ¼ gi Pðt Þ þ i Pðt0 Þ ¼ 0 dt i G G  gi 0 G

ðB:1Þ

The precursor concentrations i are therefore constant and equal to Ci(t0) while the power generated by the core remains equal to P(t0). Assuming that the SCRAM is not activated at time t ¼t0 (i.e. SSCM(t0)a5, cf. Table 2) and according to Eq. (20), the derivative

812

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

of the power generated by the core at time t ¼t0 is 6 X

6 X

d Pðt Þ Pðt0 Þ ¼  bi 0 þ gi  Ci ðt0 Þ dt G i¼1 i¼1 ¼

6 X i¼1

bi

Pðt0 Þ

G

þ

6 X i¼1

gi

bi Pðt Þ ¼ 0 G  gi 0

ðB:2Þ

The power generated by the core is therefore constant and equal to P(t0) while the precursor concentrations i remains equal to Ci(t0) and the SCRAM is not activated. According to Eq. (18), the derivative of the fuel temperature in channel i at time t¼t0, with i¼1, 2, is d R  wi  Pðt0 ÞðTc,i ðt0 ÞTe Þ T ðt0 Þ ¼ dt c,i ti ðTc,i ðt0 ÞÞ R  wi  Pðt0 ÞðTe þ R  wi  Pðt0 ÞTe Þ ¼0 ¼ ti ðTe þR  wi  Pðt0 ÞÞ

ðB:3Þ

The fuel temperatures in channel i are therefore constant and equal to Tc,i(t0) while the power generated by the core remains equal to P(t0). Finally from Eqs. (B.1) to (B.3), it is concluded that the process variables Tc,i(t) with i¼1, 2, P(t), and Ci(t) with i¼1, 2, 3, 4, 5, 6, are constant and equal to Tc,i(t0), P(t0), and Ci(t0), respectively, while the SCRAM is not activated, and independently of the other variables. References [1] US Nuclear Regulatory Commission. Reactor Safety Study (WASH-1400), NUREG-75/014. Washington, DC: US Nuclear Regulatory Commission; 1975. [2] Stamatelatos M, Apostolakis G, Dezfuli H, Everline C, Guarro S, Moieni P, et al. Probabilistic risk assessment procedures guide for NASA managers and practitioners, version 1.1. Washington, DC: Office of Safety and Mission Assurance, NASA Headquarters; 2002. [3] Siu N. Risk assessment for dynamic systems: an overview. Reliab Eng Syst Saf 1994;43(1):43–73. [4] Amendola A, Reina G. DYLAM-1: a software package for event sequence and consequence spectrum methodology, EUR 9224 EN. Luxembourg: Commission of European Communities; 1984. [5] Cacciabue PC, Amendola A, Cojazzi G. Dynamic logical analytical methodology versus fault tree: the case study of the auxiliary feedwater system of a nuclear power plant. Nucl Technol 1986;74(2):195–208. [6] Cojazzi G. The DYLAM approach for the dynamic reliability analysis of systems. Reliab Eng Syst Saf 1996;52(3):279–96. [7] Acosta C, Siu N. Dynamic event trees in accident sequence analysis: application to steam generator tube rupture. Reliab Eng Syst Saf 1993;41(2):135–54. [8] Hsueh KS, Mosleh A. The development and application of the accident dynamic simulator for dynamic probabilistic risk assessment of nuclear power plants. Reliab Eng Syst Saf 1996;52(3):297–314. [9] Kloos M, Peschke J. MCDET: a probabilistic dynamics method combining Monte Carlo simulation with the discrete dynamic event tree approach. Nucl Sci Eng 2006;153(2):137–56. [10] Hakobyan A, Aldemir T, Denning R, Dunagan S, Kunsman D, Rutt B, et al. Dynamic generation of accident progression event trees. Nucl Eng Des 2008;238(12):3457–67. [11] Pickard, Lowe and Garrick Inc. Seabrook station probabilistic safety assessment, PLG-0300. Newport Beach: (prepared for) Public Service Company of New Hampshire and Yankee Atomic Electric Company; 1983. [12] Swaminathan S, Smidts C. The event sequence diagram framework for dynamic probabilistic risk assessment. Reliab Eng Syst Saf 1999;63(1):73–90. [13] Swaminathan S, Smidts C. An application of the ESD framework to the probabilistic risk assessment of dynamic systems. In: Kondo S, Furuta K, editors. Proceedings of the fifth international conference on probabilistic safety assessment and management, Osaka, Japan, vols. 1–4. Tokyo: Universal Academy Press; 2000. p. 1283–9. [14] Matsuoka T, Kobayashi M. GO-FLOW: a new reliability analysis methodology. Nucl Sci Eng 1988;98(1):64–78. [15] Matsuoka T, Kobayashi M. An analysis of a dynamic system by the GO-FLOW methodology. In: Cacciabue C, Papazoglou IA, editors. Proceedings of the European safety and reliability conference/third international conference on probabilistic safety assessment and management, Crete, Greece. London: Springer; 1996. p. 1547–52. [16] Cepin M, Mavko B. A dynamic fault tree. Reliab Eng Syst Saf 2002;75(1): 83–91. [17] Bouissou M, Bon JL. A new formalism that combines advantages of fault-trees and Markov models: Boolean logic driven Markov processes. Reliab Eng Syst Saf 2003;82(2):149–63.

[18] Distefano S, Puliafito A. Dynamic reliability block diagrams VS dynamic fault trees. In: IEEE, editor. Proceedings of the 52nd annual reliability and maintainability symposium, Orlando, USA. New York: IEEE; 2006. p. 71–6. [19] Shin SK, Seong PH. Review of various dynamic modeling methods and development of an intuitive modeling method for dynamic systems. Nucl Eng Technol 2008;40(5):375–86. [20] Garrett CJ, Guarro SB, Apostolakis GE. The dynamic flowgraph methodology for assessing the dependability of embedded software systems. IEEE Trans Syst Man Cybern 1995;25(5):824–40. [21] Yau M, Apostolakis G, Guarro S. The use of prime implicants in dependability analysis of software controlled systems. Reliab Eng Syst Saf 1998;62(1–2): 23–32. [22] Dutuit Y, Chatelet E, Signoret JP, Thomas P. Dependability modelling and evaluation by using stochastic Petri nets: application to two test cases. Reliab Eng Syst Saf 1997;55(2):117–24. [23] Kermisch C, Labeau PE, Lardeux E, Chabot JL. Implementation of hybrid Petri nets—lessons learned from their application to a SMR unit. In: Spitzer C, Schmocker U, Dang VN, editors. Proceedings of the seventh international conference on probabilistic safety assessment and management/European safety and reliability conference, Berlin, Germany, vols. 1–6. London: Springer; 2004. p. 681–6. [24] Sknourilova P, Bris R. Coloured Petri nets and a dynamic reliability problem. Proc Inst Mech Eng Part O J Risk Reliab 2008;222(4):635–42. [25] Aldemir T. Computer-assisted Markov failure modeling of process control systems. IEEE Trans Reliab 1987;R-36(1):133–44. [26] Aldemir T. Utilization of the cell-to-cell mapping technique to construct Markov failure models for process control systems. In: Apostolakis G, editor. Proceedings of the first international conference on probabilistic safety assessment and management, Beverly Hills, USA, vol. 2. New York: Elsevier; 1991. p. 1431–6. [27] Tombuyses B, Aldemir T. Dynamic PSA of process control systems via continuous cell-to-cell mapping. In: Cacciabue C, Papazoglou IA, editors. Proceedings of the European safety and reliability conference/third international conference on probabilistic safety assessment and management. Crete, Greece. London: Springer; 1996. p. 1541–6. [28] Zhang H, Dufour F, Dutuit Y, Gonzalez K. Piecewise deterministic Markov processes and dynamic reliability. Proc Inst Mech Eng Part O J Risk Reliab 2008;222(4):545–51. [29] Lair W, Ziani R, Mercier S, Roussignol M. Modeling and quantification of aging systems for maintenance optimization. In: IEEE, editor. Proceedings of the 56th annual reliability and maintainability symposium, San Jose, USA. New York: IEEE; 2010. [30] Smidts C, Devooght J. Probabilistic reactor dynamics—II: a Monte Carlo study of a fast reactor transient. Nucl Sci Eng 1992;111(3):241–56. [31] Marseguerra M, Zio E. Monte Carlo approach to PSA for dynamic process systems. Reliab Eng Syst Saf 1996;52(3):227–41. [32] Labeau PE. Monte Carlo estimation of generalized unreliability in probabilistic dynamics.1. Application to a pressurized water reactor pressurizer. Nucl Sci Eng 1997;126(2):131–45. [33] Perez Castaneda GA, Aubry JF, Brinzei N. Performance assessment of systems including conflict in the context of dynamic reliability. Int J Adapt Innovative Syst 2010;1(3–4):233–47. ¨ [34] Aldemir T, Siu N, Cacciabue P, Goktepe B, Molesh A. Reliability and safety assessment of dynamic process systems. NATO ASI series, series F: computer and system sciences. Berlin: Springer; 1994. [35] Labeau PE, Smidts C, Swaminathan S. Dynamic reliability: towards an integrated platform for probabilistic risk assessment. Reliab Eng Syst Saf 2000;68(3):219–54. [36] Aldemir T, Stovsky MP, Kirschenbaum J, Mandelli D, Bucci P, Mangan LA, et al. Dynamic reliability modeling of digital instrumentation and control systems for nuclear reactor probabilistic risk assessments, NUREG/CR-6942. Washington, DC: US Nuclear Regulatory Commission; 2007. [37] Kirschenbaum J, Bucci P, Stovsky M, Mandelli D, Aldemir T, Yau M, et al. Benchmark system for comparing reliability modeling approaches for digital instrumentation and control systems. Nucl Technol 2009;165(1): 53–95. [38] Aldemir T, Guarro S, Mandelli D, Kirschenbaum J, Mangan LA, Bucci P, et al. Probabilistic risk assessment modeling of digital instrumentation and control systems using two dynamic methodologies. Reliab Eng Syst Saf 2010; 95(10):1011–39. [39] Devooght J, Smidts C. Probabilistic reactor dynamics–I: the theory of continuous event trees. Nucl Sci Eng 1992;111(3):229–50. [40] Davis MHA. Piecewise-deterministic Markov processes—a general class of non-diffusion stochastic models. J R Stat Soc 1984;46(3):353–88. [41] Davis MHA. Markov models and optimization. Monographs on statistics and applied probability series, vol. 49. London: Chapman and Hall; 1993. [42] Devooght J, Smidts C. Probabilistic dynamics as a tool for dynamic PSA. Reliab Eng Syst Saf 1996;52(3):185–96. [43] Devooght J, Labeau PE. Moments of the distributions in probabilistic dynamics. Ann Nucl Energy 1995;22(2):97–108. [44] Labeau PEA. Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics. Reliab Eng Syst Saf 1996;52(1):65–75. [45] Cocozza-Thivent C, Eymard R, Mercier S, Roussignol M. Characterization of the marginal distributions of Markov processes used in dynamic reliability. J Appl Math Stoch Anal 2006:1–18.

F. Brissaud et al. / Reliability Engineering and System Safety 96 (2011) 793–813

[46] Cocozza-Thivent C, Eymard R, Mercier S. A finite volume scheme for dynamic reliability models. IMA J Numer Anal 2006;26(3):446–71. [47] Devooght J, Smidts C. Probabilistic reactor dynamics—III: a framework for time-dependent interaction between operator and reactor during a transient involving human error. Nucl Sci Eng 1992;112(2):101–13. [48] Aldemir T, Miller DW, Stovsky MP, Kirschenbaum J, Bucci P, Fentiman AW, et al. Current state of reliability modeling methodologies for digital systems and their acceptance criteria for nuclear power plant assessments, NUREG/ CR-6901. Washington, DC: US Nuclear Regulatory Commission; 2006. [49] Barger P, Thiriet JM, Robert M. Dynamical reliability and availability evaluation and validation of distributed control systems. In: IEEE, editor. Proceedings of the 19th IEEE instrumentation and measurement technology conference, Anchorage, USA, vol. 1. New York: IEEE; 2002. p. 837–42. [50] Ghostine R, Thiriet JM, Aubry JF. Variable delays and message losses: influence on the reliability of a control loop. Reliab Eng Syst Saf 2011;96(1): 160–71. [51] Al-Dabbagh AW, Lu L. Reliability modeling of networked control systems using dynamic flowgraph methodology. Reliab Eng Syst Saf 2010;95(11): 1202–9. [52] Zio E, Di Maio F. Processing dynamic scenarios from a reliability analysis of a nuclear power plant digital instrumentation and control system. Ann Nucl Energy 2009;36(9):1386–99. [53] Zhu D, Mosleh A, Smidts C. A framework to integrate software behavior into dynamic probabilistic risk assessment. Reliab Eng Syst Saf 2007;92(12): 1733–55. [54] Electric Power Research Institute. Guideline for performing defense-in-depth and diversity assessments for digital I&C upgrades—applying risk-informed and deterministic methods, 1002835. Palo Alto: Electric Power Research Institute; 2004. [55] Brissaud F, Barros A, Be´renguer C, Charpentier D. Dependability issues for intelligent transmitters and reliability pattern proposal. In: Bakhtadze N, Dolgui A, editors. Proceedings of the 13th IFAC symposium on information control problems in manufacturing, Moscow, Russia, vol. 13. Oxford: Elsevier; 2009. p. 298–303. [56] Brissaud F, Barros A, Be´renguer C, Charpentier D. Reliability analysis for new technology-based transmitters. Reliab Eng Syst Saf 2011;96(2):299–313. [57] International Electrotechnical Commission. IEC 60770-3, Transmitters for use in industrial-process control systems, Part 3: methods for performance evaluation of intelligent transmitters. Geneva: International Electrotechnical Commission; 2006. [58] Smith G, Bowen M. Considerations for the utilization of smart sensors. Sensors Actuators A 1995;47(1–3):521–4. [59] Brignell JE. The future of intelligent sensors: a problem of technology or ethics? Sensors Actuators A 1996;56(1–2):11–5. [60] Meijer GCM. Concepts and focus point for intelligent sensor systems. Sensors Actuators A 1994;41(1–3):183–91. [61] CIAME, Bayart M, Conrard B, Chovin A, Robert M. Intelligent sensors and actuators. Tech Ing 2005:S7520. [in French]. [62] Taner AH, Brignell JE. Aspects of intelligent sensor reconfiguration. Sensors Actuators A 1995;47(1–3):525–9. [63] Staroswiecki M. Intelligent sensors: a functional view. IEEE Trans Ind Inf 2005;1(4):238–49. [64] Blanke M, Izadi-Zamanabadi R, Bøgh SA, Lunau CP. Fault-tolerant control systems—a holistic view. Control Eng Pract 1997;5(5):693–702. [65] Brissaud F, Barros A, Be´renguer C. Improving availability and safety of control systems by cooperation between intelligent transmitters. In: Proceedings of the 10th international conference on probabilistic safety assessment and management, Seattle, USA, June 7–11, 2010. [66] Brissaud F, Smidts C, Barros A, Be´renguer C. Dynamic reliability modeling of cooperating digital-based systems. In: Ale BJM, Papazoglou IA, Zio A, editors.

[67] [68] [69]

[70]

[71]

[72]

[73] [74] [75] [76] [77] [78] [79] [80]

[81]

[82]

[83]

[84] [85]

[86] [87]

[88]

813

Reliability, risk and safety: back to the future. Boca Raton: CRC Press, Taylor & Francis Group; 2010. Aldemir T, Siu N. Reliability and safety analysis of dynamic process systems—Guest editorial. Reliab Eng Syst Saf 1996;52(3):181–3. Ross SM. Stochastic processes. 2nd ed.. USA: John Wiley & Sons; 1996. Devooght J. Uncertainty analysis in dynamic reliability. In: Mosleh A, Bari RA, editors. Proceedings of the fourth international conference on probabilistic safety assessment and management, New York City, USA, vols. 1–4. New York: Springer; 1998. p. 2263–8. Labeau PE, Amar ZO. Monte Carlo estimation of generalized unreliability in probabilistic dynamics. 2. Handling uncertainties in parameters. Nucl Sci Eng 1997;126(2):146–57. Eymard R, Mercier S. Comparison of numerical methods for the assessment of production availability of a hybrid system. Reliab Eng Syst Saf 2008;93(1):168–77. Nelson P, Wang SW. Dynamic reliability via computational solution of generalized state-transition equations for entry-time processes. Reliab Eng Syst Saf 2007;92(9):1281–93. Stefanou G. The stochastic finite element method: past, present and future. Comput Methods Appl Mech Eng 2009;198(9–12):1031–51. Murata T. Petri nets—properties, analysis and applications. Proc IEEE 1989; 77(4):541–80. Bause F, Kritzinger PS. Stochastic Petri nets—an introduction to the theory. 2nd ed.. Germany: Vieweg Verlag; 2002. Jensen K. Coloured Petri nets: basic concepts, analysis methods and practical use. 2nd ed.. London: Springer-Verlag; 2002. David R, Alla H. Petri nets for modeling of dynamic systems: a survey. Automatica 1994;30(2):175–202. Vernez D, Buchs D, Pierrehumbert G. Perspectives in the use of coloured Petri nets for risk analysis and accident modelling. Saf Sci 2003;41(5):445–63. David R, Alla H. Discrete, continuous, and hybrid Petri nets. 2nd ed. Berlin: Springer; 2010. Yang YY, Linkens DA, Banks SP. Modelling of hybrid systems based on extended coloured Petri nets. In: Antsaklis P, Kohn W, Nerode A, Sastry S, editors. Proceedings of the third workshop on hybrid systems II, New York, USA, vol. 999. Berlin: Springer-Verlag; 1995. p. 509–28. Everdij MHC, Blom HAP, Klompstra MB. Dynamically coloured Petri nets for air traffic management safety purposes. In: Papageorgiou M, Pouliezos A, editors. Proceedings of the eighth IFAC/IFIP/IFORS symposium on transportation systems, Chania, Greece, vols. 1–3. Oxford: Pergamon Press; 1997. p. 169–74. Giua A, Usai E. High-level hybrid Petri nets: a definition. In: IEEE, editor. Proceedings of the 35th IEEE conference on decision and control, Kobe, Japan, vols. 1–4. New York: IEEE; 1996. p. 148–50. Demongodin I, Koussoulas NT. Differential Petri nets: representing continuous systems in a discrete-event world. IEEE Trans Autom Control 1998; 43(4):573–9. Petri Nets World. Available from: /http://www.informatik.uni-hamburg.de/ TGI/PetriNets/S. Commission of the European Communities. Comparative analysis of a hypothetical loss-of flow accident in an irradiated LMFBR core using different computer models for a common benchmark problem, EUR 11925 EN. Luxembourg: Commission of European Communities; 1989. CPN Tools. Available from: /http://wiki.daimi.au.dk/cpntools/cpntools.wikiS. Jensen K, Kristensen LM, Wells L. Coloured Petri nets and CPN tools for modelling and validation of concurrent systems. Int J Software Tools Technol Trans 2007;9(3):213–54. CNN.com. Star survey reaches 70 sextillion. CNN.com/Science and Space 2003, Wednesday, July 23.