Petri Nets and Pseudo-Bond Graphs for a Nuclear Reactor Primary Coolant System

In this paper we demonstrate the application of timed Petri Net methodology to the frequency of initiating faults in the coolant circulation system of a ﬁctional nuclear reactor design. A second model is presented where the Petri Net is coupled with a Bond Graph representation of reactor thermodynamics, and distributions of temperatures reached by core coolant and fuel cladding during component failures is produced.


Introduction
Half of the world's operating nuclear reactors came into service between 1970and 1985(IAEA, 2019, consequently, there are many plants still in service beyond their originally intended lifespans. When these reactors were designed, the methods available to perform their risk assessment were Fault Tree (Rasmussen, 1975) and Event Tree (Andrews and Moss, 1993) analysis. However, we make use of Petri Net modeling (Petri, 1962), see figure 1, which allows one to capture any feature found in those methodologies, while also providing an easy and flexible means to represent the maintenance, repair, and aging of components. These probabilistic methodologies implicitly requires one to make assumptions about the ongoing physical conditions in the system being consid-ered. Therefore we also construct a hybrid model incorporating the Bond Graph methodology to capture the evolution of reactor core temperatures following faults arising in the primary coolant.

System Description
A demonstration system, shown schematically in figure 2, is used in this work, whose description is typical of advanced thermal reactor designs. Coolant enters the core through an inlet header and is circulated through 450 coolant channels, each containing a fuel rod clad in zircalloy-2 (Whitmarsh, 1962). On leaving the reactor, the fluid is divided into four steam separators, which extract steam to be directed to the turbine building. After condensation, feed pumps transfer coolant back to the steam separators in liquid form, from Proceedings of the 29th European Safety and Reliability Conference. which it can flow back into the core. When a controlled shutdown is required, a valve is shut on the pipe running to the turbine. This directs coolant on an alternative path into a set of eight shutdown condensers which sit in a large tank of water serving as a heat sink. The valve is serviced on a five year cycle. When the reactor shuts down, the isolation system is intended to provide decay removal for 40 days. In the event of loss of coolant pressure, four pressure accumulators sit ready to deliver high pressure emergency coolant, before transitioning to gravity fed low pressure injection from the aforementioned tank. The Petri Net model considers the sequence of events leading to a coolant fault initiating event, which can arise from either coolant loss from pipe  Places hold tokens which are given and taken by transitions when the latter's arc weights are met. If a token is placed on one the black terminal places, the simulation ends. Once a transition's firing conditions are met, it may either fire instantaneously (gray fill) or after a duration selected from a probability distribution (white fill). Arcs have a weight of 1, unless labeled otherwise. Fig. 2. Reactor core coolant circulation system schematic or steam separators, or from a critical number of concurrent feed pump failures. The system should be able to continue providing adequate coolant to the core with only two steam separators functional. If a steam separator circuit fails, a controlled shutdown is scheduled for six months time, and if a second should fail in that time, the system is shutdown immediately. Similarly the system requires that two feed pumps be running to maintain sufficient coolant, with the third standing by in case of failure, providing the necessary pumping capacity while repairs are made.

Petri Net Modeling
The Petri Net model for core coolant fault initiating events is found figure 3. Faults occurring between the turbine isolation valve and the return of coolant to the steam separators will result in an attempt to shutdown the reactor via the shutdown condensers. As this process still relies on the circulation of coolant through the steam separators, any fault in this part of the system necessitates immediate invocation of the pressure accumulators. Thus the two outcomes that are considered are the completion of the 40 day decay heat removal period, and the occurrence of an event that mandates the initiation of the emergency injection process, faults in the shutdown condenser and emergency injection systems themselves being excluded from the model at this stage. For the parameterization of timed transitions, refer to table 1.
A total of 2.6268 × 10 9 hours were simulated across 10 5 iterations of the Petri Net, giving a mean simulated duration of 3.00 years, with a median of 2.26 years. No simulations reached the limit of 100 years. A initiating event occurred in 0.85% of the simulations, with the remainder safely completing an shutdown condenser shutdown period; these simulations respectively having mean durations of 3.28 and 2.99 years. These results give a coolant initiating fault frequency of 2.84×10 −8 per operational reactor year.

Bond Graph methodology
The Bond Graph methodology (Paynter, 1961), uses the concept of the dynamical analogy to model the transfer of power through a system comprised of varying energy domains. The elements of the system, see table 2, are connected by a network of bonds, rendered on diagrams by the symbol, , whose harpoon arrowhead denotes positive directionality. These carry an "effort" variable and a "flow" variable, the product of which gives the power moving through the bond. Effort and flow are directed at junctions and transformers, where the 0-junction splits flow and enforces equal efforts over its bonds, the 1junction splits effort and enforces equal flows, and the transformer (TF) takes in effort and flow and outputs them multiplied by a constant. Once a representation of the system in Bond Graph form is realized, the set of simultaneous differential equations implied is solved to execute a simulation, in the case of this work, that is achieved in the Modelica language (Fritzson and Engelson, 1998), with use of the BondLib package (Cellier andÀngela Nebot, 2005), and the Python interface provided by JModelica (Modelon AB., 2018). For the physical model in this work, we are concerned with the transfer of thermal energy by the motion of fluid, i.e. advection, and therefore our Bond Graphs contain fluid domain and thermal domain concepts (Karnopp, 1978). However, the most convenient notation for the latter is to take temperature as the effort concept, and the enthalpy transfer rate as the flow concept. While the methodology for executing the bond graph does not change, this breaks the rule that the product of a bond's effort and flow should have the dimensionality of power, thus the thermal model is referred to a pseudo-Bond Graph. Also, the transfer of thermal energy is directly coupled to the motion of coolant; consequently, we introduce three additional elements: the pipe resistor, the heat transfer transformer, and the dual domain capacitor, the governing formulae of which are found in equations (1) to (3). A pipe resistor is always paired with a heat transfer transformer; the mass flow through the pipe resistor is used dynamically to control the heat transfer transformer's conversion value, and thus the enthalpy flow that accompanies the given motion of coolant. The Table 1. Timed transition parameters for Petri Net seen in figures 3 and 6, given in hours where applicable. The Weibull distribution is governed by a shape pameters, β and a scale parameter, η, which for a mean time to failure , t , is given by is the Gamma function. "Delay" indicates a fix duration of length, a, "uniform" indicates a uniform probability density over in (0, u]. Data sourced from Hubble and Miller (1978); Smith (1981); IAEA (1988); GRS (1990); Eide et al. (1990); Eide and Calley (1993); Cadwallader (1998); Barringer & Associates, Inc (2010), the US Nuclear Regulatory Commission, India's Atomic Energy Regulatory Board, and expert opinions. Please note that these figures are to be taken as representative values for a fictional reactor.

Transition
Distribution The behavior of the pipe resistor is governed by, (1) (Nakayama, 2013) where l P and r P are respectively the length (10 m) and radius of the pipe (0.24 m), µ is the dynamic viscosity of the fluid (8.90×10 −4 Pa·s), and ρ in is the coolant density in the tank from which it flows. ∆p andṁ are respectively the pressure differential over the pipe and the mass flow through it, these being effort and flow variables of the element, with R f being the resulting, fluidic resistance. On the Bond Graph, pipe resistors elements are given the symbol, R P .

Equations for heat transfer transformer
The transfer of thermal energy by advection in the transformer element coupled with a pipe resistor is given by,Ḣ where c p is the specific heat capacity of the coolant (4.182×10 3 J·(kg·K) −1 ) andṁ is mass flow transferred by the pipe resulting from equation (1). The 'in' subscript indicates the inwards port of the element, and the 'out' subscript indicates the resulting values on the outward port. T andḢ are respectively, temperature and enthalpy flow with the aforementioned subscripts applied as relevant in equation (2). The symbol marking a heat transfer transformer is TF H .

Equations for dual-domain capacitor element
The fluidic variables, mass flow,ṁ and pressure, P , and their respective counterparts in the thermal domain, enthalpy flow,Ḣ, and temperature, T , are governed by equation (3). The form of equations (3a) to (3d) are unaltering: To limit the volume of the fluid in the tank, two governing equations are conditional, such that, where V and m are respectively the volume and mass of fluid in the tank, ρ is the density of the coolant, h is the level of the fluid relative the bottom of the tank, A T is the cross-sectional area of the tank in the horizontal plane (14.1 m 3 ), F t is used to calculated thermal expansion, for a cubic thermal expansion coefficient, β (2.07 K −1 ), with reference temperature, T ref (277.15 K), and reference density, ρ ref (1000 kg·m −3 ), F c is the additional pressure resulting from compress for a given bulk modulus, E, V T (42.3 m 3 ), is the volume of the tank, and g (9.81 N·kg −1 ) is the local gravitational field strength.
Internally, the two energy domains are coupled via the specific heat capacity, such that, On the Bond Graph, the tank element is given the symbol, C T .

Core temperature model
The pseudo-Bond Graph model of the thermofluid dynamics of core coolant throughput is seen in figure 4. Heat transfer between the cladding and the core coolant is governed by the total thermal resistance, R thm , where R thc , R rad , and R cnv are the contributions from thermal contact transfer, radiation, and forced convection respectively, such that, (Holman, 1990) where r 1 is the radius of outer fuel rod including cladding (0.005 m), r 2 is the inter radius of a fuel assembly pressure tube (0.1 m), A x is the rod surface area in contact with coolant (848 m 3 for 450 3 m long fuel rods), and k 1 and k 2 are respectively the thermal conductivities of the cladding (14.19 W·(m·K) −1 Whitmarsh (1962)) and coolant (0.68 W·(m·K) −1 Incropera and de Witt (1990)),  (Holman, 1990) where h c is the convection heat transfer coefficient (taken to be 5 × 10 4 W·(m 2 ·K) −1 for force convection of water in a pipe (Whitelaw, 2011)).
At full power, the reactor produces 9.2×10 8 W in heat, on shutdown its output is cut to 6%. From that point, core decay heat output falls exponentially, reaching 1.5% at one hour after shutdown. Thus enthalpy entering the Bond Graph system from nuclear processes,Ḣ N , at time, t, is given by,Ḣ whereḢ full is full reactor output, t SI is the time at which shutdown is initiated, and λ is a constant equal to ( ln (4) /3600) s −1 . The coolant heat in is the product of the mass flow into the coreṁ CIn , the specific heat capacity of the coolant, and its temperature, T CIn , such that, The heat of the coolant flowing out of the system is the temperature of coolant in the core. The outgoing pressure, P COut , is, where P A and P core are respectively atmospheric pressure (101325 Pa) and the pressure of coolant in the core. The cladding has a thermal heat capacity of 2.44×10 5 J·K −1 based on a total mass and specific heat capacity of 697 kg and 347 J·(kg·K) −1 respectively, with an emissivity of 0.433 (Whitmarsh, 1962;Murphy and Havelock, 1976).

Results in isolation
Before combining the Bond Graph with the Petri Net, two simulations were conducted to observe its behavior running alone. In the first, the initial temperatures core coolant and cladding were 303 K and the model ran until equilibrium was reached, settling at 500 K for the core coolant, and 933 K for the cladding. In the second, a 2.25 hour long execution of the Bond Graph was conducted, with the initial temperatures at equilibrium. At 0.25 hours into the simulation, core shutdown was initiated, and the extraction of decay heat was observed. The results are shown in figure 5; two hours after shutdown, the cladding and core coolant have cooled to be close to the temperature of incoming coolant. In both simulations, a constant coolant inflow rate of 2200 kg·s −1 was used throughout.

Interaction between models
To control the Bond Graph, one must consider the physical consequences of faults in the Petri Net, and what physical consequences arising in the Bond Graph constitute a failure. To this end we extend the Petri Net to produce a severity of break associated with coolant circulation faults, with the resulting information transmitted to the Bond Graph. Likewise, the Petri Net structure is also extended such that places exist to deposit tokens when the Bond Graph enters a state that requires action in the Petri Net. A diagram illustrating these additions and their relationship to the Bond Graph is found in figure 6, where one will see the reduction of coolant throughput resulting from ruptures or feed pump failure, such that the incoming coolant mass flow,ṁ CIn , is,  0.5 if only one pump is functional, and 0 is all pumps have failed, and d SS,i , d OP , and d TP are damaged factors for the steam separators, other piping (inlet header, coolant channels, and first section of pipe to turbine), and piping beyond the turbine isolation valve (the latter two sections of the turbine pipe, the turbine itself, and the condenser), respectively, and all of which represent the coolant lost to a breach, taking a value between 0 (no rupture) and 1 (complete loss of coolant circulation through component). If the cladding's temperature reaches 60% of its melting point, 2122 K (Whitmarsh, 1962), a shutdown command is given, and at 80%, the simulation terminates on the presumption that the reactor will imminently sustain damage as a result of the extreme temperature.

Hybrid model results
5000 executions of the joint Petri Net-Bond Graph model were conducted. When simulations encountered a problem, the vast majority were able to retain control the core temperature by going to shutdown. However, in some cases, the coolant fault that emerged was too severe and the simulation ended with the core coolant and cladding at elevated temperatures, with one instance of termination resulting from the cladding temperature surpassing the 60% of melting temperature threshold. In figure 7, a histogram is plotted show- ing the distribution of core coolant and cladding temperatures, both at the end of a simulation and the highest intermediate values reached. When temperatures rise above the normal operating values, the increase is generally small, with the highest core coolant temperature of 4402 simulations falling into the 505 to 555 K range, and the highest cladding temperature of 4572 simulations being within 935 to 975 K.
As few simulations ended in a fault, a second set of 5000 simulations was conducted. In this set, the ability to close the isolate valve and shutdown the reactor was removed, and therefore these results can be considered a deeper exploration of the subset of systems encountering such a fault. The distribution of final and highest intermediate temperatures reached for these results is shown in figure 7. Unsurprising, more simulations are seen to get significantly hot when the reactor shutdown procedure is eliminated, although only three instances are seen in which the cladding temperature passes the 60% threshold. In exactly 1000 simulations, the highest temperature reached by both the core coolant and cladding was 50 K above prefault levels, with 440 simulations experiencing increases of 100 K.

Discussion
The Petri Net constructed in this work has given an estimation of the frequency of faults requiring activation of emergency core coolant injection. A subset of which of these would experience further faults, leading to a loss of coolant event. Our continuing work will model the sequence that follows an initiating event, using the data produced in this work to calculated overall failure frequency, and consideration of each the major reactor subsystems will be given before a whole system depiction is constructed. A joint Petri Net-Bond Graph model was presented for the demonstration system, in which one observes the range of temperatures reached by core coolant and cladding during the occurrence of faults. In these simulations, elevated cladding temperatures were seen, but no instances in which melting would be imminent arose. Where applicable, we will continue to couple our Petri Nets with a Bond Graph representation of system thermodynamics as the scope of our studies expands. It is apparent that a greater volume of simulations need be conducted to capture the full spectrum of outcomes.