Abstract

This paper presents a model of a simplified boiling water reactor (SBWR) to analyze the steady-state and transient behavior. The SBWR model is based on approximations of lumped and distributed parameters to consider neutronics and natural circulation processes. The main components of the model are vessel dome, downcomer, lower plenum, core (channel and fuel), upper plenum, pressure, and level controls. Further consideration of the model is the natural circulation path in the internal circuit of the reactor, which governs the safety performance of the SBWR. To demonstrate the applicability of the model, the predictions were compared with plant data, manufacturer_s predictions, and RELAP5 under steady-state and transient conditions of a typical BWR. In steady-state conditions, the profiles of the main variables of the SBWR core such as superficial velocity, void fraction, temperatures, and convective heat transfer coefficient are presented and analyzed. The transient behavior of SBWR was analyzed during the closure of all main steam line isolation valves (MSIVs). Our results in this transient show that the cooling system due to natural circulation in the SBWR is around 70% of the rated core flow. According to the results shown here, one of the main conclusions of this work is that the simplified model could be very helpful in the licensing process.

1. Introduction

Recently, a joint project among three Agencies: the International Energy Agency (IEA), the OECD Nuclear Energy Agency (NEA), and the International Atomic Energy Agency (IAEA) prepared the report “Innovative nuclear reactor development, opportunities for International cooperation” [1], whose main objective is the identification of opportunities to establish joint projects in developing nuclear fission reactor technologies. This study remarks that the following areas are good candidates for developing new, broad based collaborative efforts:

(i)natural circulation;(ii)high-temperature materials;(iii)passive (safety) devices;(iv)in-service inspection and maintenance methods;(v)advanced monitoring and control technologies;(vi)delivery and construction methods;(vii)safeguard technologies and approaches.

This paper focuses on the attention in natural circulation which is one of the good candidates for developing new reactor technologies. Natural convection circulation of coolant in the reactor cooling system is used to some extent in all reactor concepts, and many of these uses are innovative [2]. What is now needed is to develop or confirm the design basis and to develop and qualify computer codes to enable reliable analysis.

All the ordinary water-cooled and heavy water-cooled reactor designs rely on natural convection of the coolant to remove decay heat from the fuel after shutdown, even today’s operating water reactors also rely on natural convection to remove decay heat when forced circulation is lost.

Among the new design with natural circulation, we can mention the CAREM-25 [3] and NHR-200 [4], both use natural convection circulation of the ordinary water coolant at all power levels to remove heat from the reactor core. The energy amplifier and BREST 300 [5] use natural circulation of the lead coolant to remove decay heat from the core after shutdown, requiring the design of a reliable system to accommodate the potential solidification of the lead in regions outside the core. The energy amplifier also uses natural circulation of the coolant to remove heat produced by fission from the reactor core at all power levels; it is important to remark that this has not been done before. The FUJI [2] uses natural circulation of the molten salt reactor coolant to remove decay heat from the core after shutdown. Therefore, the ability to predict the performance of systems operating in natural convection circulation mode has fundamental importance.

One of the main design concepts is presently proposed by general electric (GE): the simplified boiling water reactor (SBWR) [6] and the economic and simplified boiling water reactor (ESBWR) [7]. According to the manufacturer, these new design reactors working with natural circulation provide simplification over previous BWRs. The improvements are accomplished by the removal of the recirculation pumps and associated motors, piping, valves, heat exchanges, control, and electrical support systems that exist with forced recirculation, eliminating flow disturbance by abnormal pump behavior. The main differences between natural and forced circulations are the additions of an exhaust pipe above the reactor core to stabilize and direct the steam and water flow above the core, and the correspondingly open and downcomer annulus that reduces the flow resistance and provides additional driving head pushing the water to the bottom of the core. Among the benefits of these designs are simplification of the operation, higher plant availability, 20% reduction in the cost of operation and maintenance, and the reduction of personnel dose by the elimination of recirculation pumps maintenance.

Simplified and convenient models are used in computer codes to perform analysis in early stages of the design of experimental reactor facilities and Generation IV reactors. The work of Tian et al. [8] showed a simplified neutronic and thermalhydraulic model used in the China Advanced Research Reactor (CARR) and its application with a simulation of loss of power accident without emergency cooling system. Other application is the engineering safety analysis for regulatory licensing process of light water reactors (LWR) that were performed in the past using simplified one-dimensional models of thermalhydraulics systems with the appropriate initial and boundary conditions for transient analysis in two-phase flow phenomena due that it can provide overestimated results.

This paper presents an SBWR model for transient analysis, where the thermalhydraulic model that describes the dynamic behavior of the lower and upper plena, and the reactor core, as well as fuel temperature model were based on distributed parameter approximations. The vessel dome, downcomers, recirculation loops, and neutron process models were based on lumped parameter approximations.

A multinode fuel pin model is developed to describe the heat transfer process. Three regions were identified as indispensable to be considered in the heat transfer analysis: the first region corresponds to heat transfer in the fuel; the second region corresponds to heat transfer in the gap; and the third region corresponds to heat transfer in the clad, whose temperatures are determined by the rate of heat convection due to core flow.

Reactor power is calculated from a point reactor kinetics model with six groups of delayed neutrons. The reactivity due to Doppler effect, void fraction, moderator temperature, and control rod drive was considered in this model.

This paper is organized as follows. The SBWR system description is presented in Section 2; mathematical models are described in Section 3; and Section 4 shows the results and discussions. The qualification assessment is presented in the appendix.

2. System Description

The SBWR configuration and the flow paths are illustrated in Figure 1. The reactor natural circulation loop circulates the required coolant flow through the reactor core. The flow within the reactor vessel provides a continuous internal natural circulation path for a major portion of the core coolant flow. The core flow is taken from the vessel and discharged into the lower core plenum. The coolant water passes along the individual fuel rods inside the fuel channel where it boils and becomes a two-phase steam/water mixture. In the core, the two-phase fluid generates upward flow through the axial steam separators, while the steam continues through the dryers and flows directly out through the steam lines into the turbine generator. The condensate flow is then returned through the feed water heaters by the condensate feed water pumps into the vessel. Finally, the water, which is separated from the steam in the steam separators, flows downward in the periphery of the reactor vessel and mixes with the incoming main feed flow from the turbine.

Table 1 shows the nominal values used in the simulation, which correspond to design characteristics of the SBWR [9].

3. Mathematical Model

Figure 2 is a schematic diagram of the SBWR where the arrangement of the computational cells of the SBWR model is shown. The reactor vessel was divided into five zones. Two of these zones, the vessel dome and the downcomer, have variable volume according to the vessel water level. The three fixed volume zones are the lower plenum, the upper plenum and steam separators, and the reactor core. Due to its importance on the model performance, the latter was subdivided into twelve one-dimensional nodes. The reactor model is completed by including natural recirculation loop as coolant, neutron kinetics, fuel rod temperature, and control models. In addition, this model uses a set of empirical correlations valid for the normal range of SBWR operating conditions. The model derivation is based in lumped and distributed parameters approximation; Table 2 presents the submodels and approximations used.

The following assumptions were adopted in the development of the SBWR model:

(i)the point-kinetic equations are used to describe the evolution of neutron population and precursor concentration;(ii)delayed neutrons are modeled using six groups;(iii)fuel rod is described by a distributed parameters model, which includes fuel, gap, and cladding. This model is crucial due that velocity of heat transfer from fuel to the coolant has important impact during transient behavior;(iv)the multiple parallel channels in SBWR are lumped into one average channel. The drift-flux model is used to consider velocity difference between liquid and vapor phases;(v)the axial power profile is modeled as a superposition of a profile space-dependent and a space-independent term, that is, time-dependent;(vi)the downcomer section is assumed to be filled with single-phase fluid. The thermohydraulics system itself consists of five different sections: lower plenum, core, upper plenum, stem separators, and downcomer. The heated section is the core, and the others are unheated sections;(vii)the gas phase is saturated at the vessel pressure.

3.1. Thermalhydraulic Model

The thermalhydraulic model consists of five balance equations: (1) liquid-and gas-phase mass balances, (2) mixture momentum, (3) mixture energy, (4) liquid-phase energy, together with a drift flux approach [10], for the analysis of phase separation. The nonequilibrium two-phase flow for the volumetric vapor generation rate in subcooled boiling was considered using Saha and Zuber approximation [11, 12].

Nonequilibrium two-phase flow conditions are present in all two-phase flows, and the effects of nonequilibrium states can be important at steady-state and off-normal conditions encountered in operational transients. In boiling water reactor, for example, vapor voids are formed in the core when the bulk liquid is below the saturation state; this phenomena is known as subcooled boiling. These voids have a very strong effect relative to the power distribution in the core both at quasi-steady conditions and during operational transient.

The thermalhydraulic model is used to describe the dynamic behavior of the lower and upper plena, downcomer, and the reactor core. The five-equation model was based on the following assumptions: (i) one-dimensional flow, (ii) no cross flow, and (iii) kinetic and potential energy effects are neglected in the energy equation.

The model is based on liquid- (denoted by subscript 𝑙) and gas- (denoted by subscript 𝑔) phase mass and energy balances, and a mixture momentum balance (denoted by subscript 𝑚): 𝜕𝜌𝜕𝑡𝑔𝜀𝑔+𝜕𝑗𝜕𝑧𝑔𝜌𝑔𝜕=Γ,𝜌𝜕𝑡𝑙𝜀𝑙+𝜕𝑗𝜕𝑧𝑙𝜌𝑙=Γ,(1)𝜕𝜌𝜕𝑡𝑙𝑙𝜀𝑙𝜀𝑙𝜕𝑝+𝜕𝜕𝑡𝜌𝜕𝑧𝑙𝑙𝑗𝑙=𝑞𝑙𝑃𝐻𝐴𝑥𝑠+𝑞𝑙𝜀𝑙Γ𝑓,(2)𝜕𝜌𝜕𝑡𝑚𝑚𝜕𝑝𝜕𝑡+𝐺𝑚𝜕𝑚=𝑞𝜕𝑧𝑃𝐻𝐴𝑥𝑠+𝑞,(3)𝜕𝐺𝑚+𝜕𝜕𝑡𝐺𝜕𝑧2𝑚𝜌𝑚=𝜕𝑝𝜕𝑧𝑓𝐺2𝑚2𝜌𝑚𝐷𝑔𝜌𝑚,(4) where 𝜌𝑚=𝜌𝑔𝜀𝑔+𝜌𝑙𝜀𝑙,𝜌(5)𝑚𝑚=𝜌𝑔𝑔𝜀𝑔+𝜌𝑙𝑙𝜀𝑙,𝐺(6)𝑚=𝜌𝑔𝑗𝑔+𝜌𝑙𝑗𝑙,𝜀(7)𝑔+𝜀𝑙=1.(8)

In these equations, 𝜌 is density, 𝜀𝑔 is the void fraction; 𝑗 is superficial velocity; Γ is the volumetric vapor generation rate; is the enthalpy; 𝑝 is the pressure; 𝑞 is the heat flux; 𝑃𝐻 is the heated perimeter; 𝐴𝑥𝑠 is the cross-sectional area; 𝑓 is the friction factor; 𝐷 is the hydraulic diameter; and 𝑞 is the volumetric heat generation by gamma radiation. For unheated sections (e.g., downcomer, lower and upper plenums) 𝑞=0.

The equations of state are usually expressed in terms of two independent variables in the form 𝜌𝑙=𝜌𝑙(𝑝,𝑙) for density and 𝑇𝑙=𝑇𝑙(𝑝,𝑙) for temperature. If the fluid is at the saturation state, the state properties are functions of the pressure only, and the equations of state can be written as 𝜌𝑘=𝜌𝑘(𝑝) to density and 𝑇𝑘=𝑇𝑘(𝑝) where 𝑘=𝑙 or 𝑔. For the gas phase and in all cases 𝑔=𝑔(𝑝).

Depending on the mixture enthalpy 𝑚, the following two cases were considered.

Case 1. If 𝑓(𝑝)<𝑚<𝑔(𝑝), the gas and liquid phases are saturated at reactor pressure and 𝑙=𝑓(𝑝).

Case 2. If 𝑚<𝑓 and 𝑙<𝑓, liquid phase occurs first and, subsequently, subcooled boiling appears.

Here, 𝑓 and 𝑔 are the enthalpy of liquid and vapor saturated at the vessel pressure, respectively.

For bulk boiling condition, both the gas and liquid phases are in saturated condition (thermodynamic equilibrium condition) at the vessel pressure, then if 𝑙=𝑓(𝑝); 𝜌𝑙=𝜌𝑓(𝑝) in (3) and combining it with (1), the volumetric vapor generation rate Γ is given by: 1Γ=𝑓𝑔[𝑞𝑃𝐻𝐴𝑥𝑠+𝑞+1𝜌𝑔𝜀𝑔𝜕𝑔𝜕𝑝𝜌𝑓𝜀𝑙𝜕𝑓𝜕𝑝𝜕𝑝𝜕𝑡𝜌𝑓𝑗𝑙𝜕𝑓𝜕𝑧].(9)

In these equations, 𝜕𝑔/𝜕𝑧=0, due that the gas phase is saturated at the vessel pressure, 𝑓𝑔 is the difference in specific enthalpy of saturated vapor and liquid. For subcooled boiling, the volumetric vapor generation is given by [11]: 1Γ=𝑓𝑔[𝑞𝑃𝐻𝐴𝑥𝑠1𝑓𝑙𝑓𝑙𝑑+𝑞𝐻𝑜𝑓𝑔𝜈𝑓𝑔𝜀𝑔𝑇𝑠𝑇𝑙+𝜀𝑔1𝜌𝑔𝜕𝑔𝜕𝑝𝜕𝑝𝜕𝑡],(10)where 𝐻𝑜 is the condensation parameter, and the liquid subcooling (𝑓𝑙𝑑) at the initiation of subcooled boiling is given by [12]𝑓𝑙𝑑=0.0022𝑞𝐷𝐶𝑝𝑙/𝑘𝑙,Pe<70,000,154𝑞/𝐺,Pe>70,000.(11)

Here, 𝑘𝑙 is the liquid thermal conductivity; 𝐶𝑝𝑙 is the liquid specific heat; and 𝑇𝑠 is the saturation temperature.

The mass flow rate of each phase is, respectively, given by 𝑊𝑔=𝑗𝑔𝐴𝑥𝑠𝜌𝑔,𝑊𝑙=𝑗𝑙𝐴𝑥𝑠𝜌𝑙,(12)where 𝑗𝑙=𝑗𝑗𝑔. In order to calculate the amount of steam and liquid from mixture properties, the drift-flux model was used [10]. This model provides an algebraic relationship between the relative velocities of the flowing liquid and steam and the void fraction. The basic equation of the drift flux model is as follows: 𝑗𝑔=𝐶𝑜𝑗+𝜈𝑔𝑗𝜀𝑔.(13)

The distribution parameter 𝐶𝑜 and the average drift velocity 𝜈𝑔𝑗 are represented by a correlation for each of the two-phase flow regimes, bubbly and slug in the present case, in vertical pipes. The superficial velocity is calculated by: 1𝑗=[𝜌𝑔1𝜌𝑙𝜀Γ𝑔𝜌𝑔𝜕𝜌𝑔+𝜀𝜕𝑝𝑙𝜌𝑙𝜕𝜌𝑙𝜕𝑝𝜕𝑝𝜕𝑡𝜉sub+𝑊]Δ𝑧gin𝜌𝑔+𝑊lin𝜌𝑙1𝐴𝑥𝑠,(14)where𝜉sub=𝜀𝑙𝜌𝑙𝜕𝜌𝑙𝜕𝑙𝜕𝑙𝜕𝑡,Subcooledboiling,0,Thermodynamicequilibrium.(15)Here, Δ𝑧 is the cell length in axial direction; 𝑊gin and 𝑊lin are the inlet mass flow rate of the gas and liquid phases, respectively.

Equations (1)–(15) represent the thermalhydraulic model, which consider nonequilibrium two-phase flow effects and velocities difference between liquid and vapor phases. The main advantage of this model is its simplicity requiring only and few constitutive equations. In fact, only well-known, wall-to-mixture heat transfer and friction constitutive relationships are required. The term 𝜕𝑝/𝜕𝑡 in the equations of the thermalhydraulic model corresponds to vessel dome, which is presented in Section 3.5 The thermalhydraulic parameters are presented in Table 3.

For time discretization, an explicit Euler method was used, with a fixed time step of 0.01 second to solve the equations set of the thermalhydraulic models. This time step is highly recommended for transient with fast change in vessel pressure. Regarding the spatial discretization, the scheme used is known as a first-order, upwind finite difference scheme with a truncation error of order Δ𝑧.

3.2. Reactor Power Model

The reactor power is given by 𝑃(𝑡,𝑧)=𝑛(𝑡)𝐹(𝑧)𝑃0,(16)where 𝐹(𝑧) is the axial power factor (which is given by Figure 3); 𝑃0 is nominal power; and 𝑛(𝑡) is the normalized neutron flux, which is calculated from a point reactor kinetics model with six groups of delayed neutrons: 𝑑𝑛(𝑡)=𝑑𝑡𝜌(𝑡)𝛽Λ𝑛(𝑡)+6𝑖=1𝜆𝑖𝑐𝑖(𝑡),(17)𝑑𝑐𝑖(𝑡)=𝛽𝑑𝑡𝑖Λ𝑛(𝑡)+𝜆𝑖𝑐𝑖(𝑡),𝑖=1,2,,6,(18)where 𝑐𝑖 is a delayed neutron concentration of the 𝑖th precursor group normalized with the steady-state neutron density; 𝜌 is the net reactivity; 𝛽 is the neutron delay fraction; Λ is the neutron generation time; and 𝛽𝑖 is the portion of neutrons generated by the 𝑖th group. The initial conditions are given by 𝑛(0)=𝑛0 and 𝑐𝑖(0)=𝛽𝑖𝑛0/Λ𝜆𝑖. The parameters of the kinetics model are presented in Table 4.

The net reactivity of the nuclear reactor includes four main components: feedback reactivity due to the void fraction in the two-phase flow 𝜌𝑣, Doppler effect 𝜌𝐷 due to fuel temperature, moderator temperature 𝜌𝑚, and reactor control rods 𝜌cr. Therefore, we express the feedback reactivities as: 𝜌(𝑡)=𝜌𝑣(𝑡)+𝜌𝐷(𝑡)+𝜌𝑚(𝑡)+𝜌cr(𝑡).(19)

The kinetics point equations are stiff in the coefficients because they differ in several orders of magnitude, especially in (17). The variable implicit integration method [14] was used to solve (17), and the Euler method in an explicit form was used to solve the delayed precursor concentration given by (18). The reactivity correlations are presented in Table 5.

3.3. Fuel Model

The fuel assembly temperature distribution was obtained considering each radial node at each of the twelve hydraulic axial nodes in the core. Two radial nodes are considered for the clad and the gap, two more nodes for the boundary condition evaluation, and four nodes are considered for one equivalent fuel element, as illustrated in Figure 4.

A detailed multinode fuel pin model is developed for this study. The fuel heat transfer formulation is based on the following fundamental assumptions: (i) axis-symmetric radial heat transfer; (ii) the heat conduction in the axial direction is negligible in respect to the heat conduction in the radial direction; (iii) the volumetric heat rate generation in the fuel is uniform in each radial node; and (iv) storage of heat in the fuel cladding and gap is negligible.

Under these assumptions, the transient temperature distribution in the fuel pin, initial, and boundary conditions are given by: (𝜌𝐶𝑝)𝜕𝑇1𝜕𝑡=𝑘𝑟𝜕𝑟𝜕𝑟𝜕𝑇𝜕𝑟+𝑞(𝑡),at𝑟0𝑟𝑟6,(20)I.C.𝑇(𝑟,0)=𝑇(𝑟),at𝑡=0,(21)B.C.1𝑘𝜕𝑇𝜕𝑟=𝐻𝑇𝑇𝑚,at𝑟=𝑟6,(22)B.C.2𝜕𝑇𝜕𝑟=0,at𝑟=𝑟0.(23)

In these equations, 𝑟 is the cylindrical radial coordinate; 𝑟0 and 𝑟6 are defined in Figure 4, 𝑞(𝑡)=𝑃(𝑡)/𝑉𝑓 at each axial node, where 𝑃 is given by (16); 𝑇𝑚 is the moderator temperature; and 𝐻 is the convective heat transfer coefficient, which is presented in Section 3.4.

In the present solutions, the fuel element is represented by a one-dimensional, mesh-centred grid consisting of a variable number of radial elements at each axial position. The differential equations described previously are transformed into discrete equations using the control volume formulation technique in an implicit form [15].

Application of the control volume formulation enables the equations for each region (fuel, gap, and clad) to be written as a single set of algebraic equations for the sweep in the radial direction: 𝑎𝑗𝑇𝑡+Δ𝑡𝑗1+𝑏𝑗𝑇𝑗𝑡+Δ𝑡+𝑐𝑗𝑇𝑡+Δ𝑡𝑖,𝑗+1=𝑑𝑗,(24) where 𝑇𝑡+Δ𝑡𝑗1, 𝑇𝑗𝑡+Δ𝑡, and 𝑇𝑡+Δ𝑡𝑗+1 are unknowns; 𝑎𝑗, 𝑏𝑗, 𝑐𝑗, and 𝑑𝑗 are coefficients, which are computed at the time 𝑡. When these equations are put into a matrix form, the coefficient matrix is tridiagonal. The solution procedure for the tridiagonal system is the Thomas algorithm [15], which is the most efficient algorithm for this type of matrices. The coefficients 𝑎𝑗, 𝑏𝑗, 𝑐𝑗 are dependents of thermophysical properties, that is, thermal conductivity, density, and specific heat, and since they are at function of 𝑇𝑗𝑡+Δ𝑡, at least one iteration is needed.

The interaction of the model with the neutron kinetics model is through evaluation of the average fuel temperature 𝑇𝑓 and the volumetric heat rate generation 𝑞. The interaction with the nuclear thermalhydraulics is through evaluation of the coolant temperature and heat transfer coefficient. The average fuel temperature is given by: 𝑇𝑓=1𝑉𝑓𝑉𝑓𝑇𝑓1𝑑𝑉𝑉𝑓𝑛𝑖=1𝑉2𝑇𝑖,2+𝑉3𝑇𝑖,3+𝑉4𝑇𝑖,4+𝑉5𝑇𝑖,5,(25) where 𝑉𝑓=𝑉2+𝑉3+𝑉4+𝑉5 is the fuel volume. Here 2,3,4, and 5 represent the radial node according to Figure 4, while that 𝑖 represents the axial index of the node. Fuel parameters are shown in Figure 4 and Table 6.

3.4. Convective Heat Transfer Coefficients

The convective heat transfer coefficient for the single-phase liquid and vapor (𝐻1𝜙) is obtained from the classical Dittus-Boelter correlation [16]. The heat transfer coefficient in both subcooled boiling and nucleate boiling regimes (𝐻2𝜙) is calculated with the Chen correlation [17]. The convective heat transfer coefficient (𝐻) is given by the maximum between 𝐻1𝜙 and 𝐻2𝜙 to avoid discontinuities, that is, 𝐻=MAX(𝐻1𝜙,𝐻2𝜙).

The interaction of the convective heat transfer coefficients with the fuel heat transfer model is through the evaluation of the clad temperature 𝑇𝑟=𝑟6 (Figure 4). The convective heat transfer coefficient 𝐻 is used in the boundary condition given by (22).

3.5. Vessel Dome Model

The vessel dome is modeled as a two-region volume, one region being liquid and the other vapor. The two regions are assumed to be at the same pressure but not necessarily at the same temperature. The dynamic model used to obtain the pressure in the vessel dome is based on macroscopic balances of mass and energy: 𝑑𝑝𝜈𝑑𝑡=𝑙𝑗𝑊𝑗𝑙+𝜕𝜈𝑙/𝜕𝑙𝑝𝑗𝑊𝑗𝑙𝑗𝑙𝑙𝑄+𝑄+𝜈𝑣𝑗𝑊𝑗𝑣+𝜕𝜈𝑣/𝜕𝑣𝑝𝑗𝑊𝑗𝑣𝑗𝑣𝑣𝑄+𝑄,(26) where 𝑄 denotes 𝑚𝑙[𝜈𝑙(𝜕𝜈𝑙/𝜕𝑙)𝑝+(𝜕𝜈𝑙/𝜕𝑝)𝑙] and 𝑄 denotes 𝑚𝑣[𝜈𝑣(𝜕𝜈𝑣/𝜕𝑣)𝑝+(𝜕𝜈𝑣/𝜕𝑝)𝑣], and 𝑊𝑗𝑙 and 𝑊𝑗𝑣 are the mass flow rates into or out of the liquid and vapor regions, respectively, and 𝑗𝑙 and 𝑗𝑣 are the enthalpy of the fluid entering or leaving the liquid and vapor regions, respectively. The mass of liquid 𝑚𝑙, mass of vapor 𝑚𝑣, enthalpy of liquid 𝑙, and enthalpy of vapor 𝑣 are given by the following balance equations: 𝑑𝑚𝑙=𝑑𝑡𝑗𝑙𝑊𝑗𝑙,𝑑𝑚𝑣=𝑑𝑡𝑗𝑣𝑊𝑗𝑣,(27)𝑑𝑙=1𝑑𝑡𝑚𝑙𝑗𝑙𝑊𝑗𝑙𝑗𝑙𝑙+𝜈𝑙𝑑𝑝,𝑑𝑡𝑑𝑣=1𝑑𝑡𝑚𝑣𝑗𝑣𝑊𝑗𝑣𝑗𝑣𝑣+𝜈𝑣𝑑𝑝.𝑑𝑡(28)The flashing and condensation flows are given by: 𝑊𝑓=𝑚𝑙(Δ𝑡𝑙𝑓𝑔𝑓𝑊),co=𝑚𝑣(Δ𝑡𝑔𝑣𝑔𝑓).(29)Here Δ𝑡 is the time step. The density and temperature in the downcomer is given by: 𝜌dw=𝜌dw𝑝,𝑙𝑇,(30)dw=𝑇dw𝑝,𝑙.(31)

The level in the vessel is calculated as a function of liquid volume, that is, 𝑁𝑙=𝑁𝑙(𝑉𝑙), where the liquid volume is given by 𝑉𝑙=𝑚𝑙𝜈𝑙.

An Euler method in an explicit form with a time step of 0.01 second was used to solve the ordinary differential equations of the vessel dome model.

3.6. Natural Circulation Loop

A natural circulation loop of the SBWR was considered in this analysis, as shown in Figure 1. Balancing the gravity head available and total loop pressure drop obtained by integration of the momentum balance leads to the model for the circulation natural. The circulation natural model includes the pressure drops and flows from the downcomer, lower and upper plenum, reactor core, and steam separators, in order to obtain the following momentum balance: 𝑑𝑊𝑐𝑑𝑡=(𝑛𝑖=1𝑙𝑖𝐴𝑖)1(𝐾psn𝑊2𝑐𝜌dw𝐾sep𝑊2sep𝜌sepΔ𝑝𝑐+Δ𝑝𝑔),(32) where (𝑙/𝐴) is the inertial term; 𝜌dw is the downcomer density (30); 𝜌sep is the steam-separator density; 𝑊sep is the flow mass through the steam separator; 𝐾psn is the support core plate loss coefficient; 𝐾sep is the separator loss coefficient; Δ𝑝𝑐 is the core pressure drop; and Δ𝑝𝑔 is the pressure drop due to gravity.

The total core pressure drop is the sum of the frictional, acceleration, and gravitational components: Δ𝑝𝑐=12𝑗=1[𝜙2𝑓0,𝑗2𝐶𝑓0𝑊2𝑐Δ𝑧𝐴2𝑥𝑠𝐷𝜌𝑓+𝑊2𝑐𝐴2𝑥𝑠1𝜌𝑐,𝑗+11𝜌𝑐,𝑗+𝜌𝑐,𝑗𝑔Δ𝑧],(33)where 𝑗=1,2,,12 are used to indicate the core node number; 𝐶𝑓0 is the single-phase friction factor; Δ𝑧𝑗 is the core node length; 𝐷 is the hydraulic diameter; 𝜌𝑓 is the density of saturated liquid; 𝜌𝑐 is the core density in each node, which is given by (5); and 𝜙2𝑓0,𝑗 is the two-phase multiplier given by [18]: 𝜙2𝑓0𝜌=1+𝑥𝑓𝜌𝑔1,(34)where 𝑥 is the vapor quality and 𝜌𝑔 is the density of the saturated vapor.

The pressure drop due to gravity without core effects, which are considered in (32), is given by: Δ𝑝𝑔=𝑔𝑁𝑙𝑉𝑙𝜌dw𝐿𝑔ps+𝐿sep𝜌sep,(35) where 𝑁𝑙, 𝐿ps, and 𝐿sep are the liquid level in the vessel, the length in the upper plenum and the length in steam separators. The parameters of this model are presented in Table 7.

The Euler method in an explicit form with a time step of 0.01 second was used to solve the ordinary differential equations (ODE) of the recirculation system model.

3.7. Control Models to Model the Control Systems of SBWR

The control models obtained considering the feed water system (FW) and the main steam line (MSL) are dummy models or auxiliary models. The interface with the process at the other end of the control loop is made by the final control element. In a vast majority of nuclear engineering processes, the final control element is an automatic control valve which throttles the flow of a manipulated variable [19].

The MSL mass flow can be calculated as [19]: 𝑊ms=𝐶𝑣ms𝑓𝐴msΔ𝑝ms𝜌𝑣|||Δ𝑝ms|||,(36)where 𝐶𝑣ms=0.4158 is the control valve size coefficient, whose value includes the MSL system; 𝑓(𝐴ms) is the fraction of the total flow area; 𝐴ms is the fraction of wide open; Δ𝑝ms=𝑝Rx6.7203×106 is the MSL pressure drop; 𝑝Rx is the reactor pressure (see Section 3.5); and 6.7203×106 is an estimated pressure. This equation allows reverse flow.

An analogous formulation is applied to calculate the FW mass flow: 𝑊fw=𝐶𝑣fw𝑔𝐴fw𝜌𝑙Δ𝑝fw,(37)where 𝐶𝑣fw=0.04595 is the control valve size coefficient, whose value includes the FW system; Δ𝑝fw=7.9976×106Pa𝑝Rx is the feed water pressure drop; and 7.9976×106Pa is an estimated pressure. The FW enthalpy is calculated as a function of the MSL mass flow rate, which was estimated using a typical BWR [20]: aa𝑊=187903.0+50281.0ms𝑊983.3ms2𝑊+8.88ms3𝑊0.03ms4,(38)where 𝑊ms=(𝑊ms/𝑊nom)×100, 𝑊ms is given by (43), and 𝑊nom is 100% of mass flow rate. The dynamic effects are considered as a response of first-order system [19]: 𝑑fw=𝑑𝑡aafw𝜏fw𝑠,(39) where 𝜏fw𝑠 is the FW time constant. Table 8 presents the time constant.

The dynamic of the valves was modeled as a response of a first-order system: 𝑑𝐴ms=𝑆𝑑𝑡cp𝐴ms𝜏ms,𝑑𝐴fw=𝑆𝑑𝑡cl𝐴fw𝜏fw,(40)where 𝜏ms and 𝜏fw are time constants of the MS valve and FW valve, respectively; 𝑆cp and 𝑆cl are control signals. The time constant values are presented in Table 8.

The proportional-integral-derivative (PID) control is used for pressure control in the reactor vessel. The proportional action changes its output signal 𝑆cp in direct proportion to the error signal 𝑒𝑝(𝑡) [19]: 𝑆cp=𝐾po+𝐾𝑝𝑒𝑝(𝑡)+𝐾𝐷𝑑𝑒𝑙𝑑𝑡+𝐾𝐼𝑒𝑝(𝑡)𝑑𝑡,(41) where 𝐾po is the polarization and the controller value output when there is no error; 𝐾𝑝, 𝐾𝐷, and 𝐾𝐼 are, respectively, the proportional, derivative, and integral control gains. The tracking error is given by: 𝑒𝑝𝑝(𝑡)=Rx7.07×106Pa𝑝nom.(42) Here, 𝑝Rx is the reactor pressure; 7.07×106Pa is the set point of pressure; and 𝑝nom is the nominal pressure. The proportional-derivative feedback control action 𝐾𝑝𝑒𝑝+𝐾𝐷(𝑑𝑒𝑙/𝑑𝑡) is used to induce certain stability margin in the controlled pressure. However, the stability of pressure does not ensure that the tracking error 𝑒𝑝(𝑡) can be made as small as desired. To this end, the integral action 𝐾𝐼𝑒𝑙(𝑡)𝑑𝑡, which contains the memory of the tracking error, is introduced to reduce the size of the steady-state tracking error.

The proportional-integral (PI) control is used for level control in the reactor vessel, which is given by [19]: 𝑆𝑐𝑙(𝑡)=𝐾𝑝𝑙𝑒𝑙(𝑡)+𝐾𝐼𝑙𝑒𝑙(𝑡)𝑑𝑡,(43) where 𝐾𝑝𝑙 is the proportional gain; 𝐾𝐼𝑙 is the integral gain; and the tracking error is given by: 𝑒𝑙𝑊(𝑡)=ms𝑊fw𝑊𝑜+𝑁inst0.97m𝑁𝑜,(44) where 𝑊𝑜 is the mass flow rate; 𝑁𝑜 is the normal level; 𝑁inst is the instrument level; and 0.97 m is the set point of the level. The proportional action moves the control valve in direct proportion to the magnitude of the error, while the integral action moves the control valve based on the time integral of the error. If there is no error, the controller output does not move. As the error goes positive or negative, the integral of the error drives the controller output either up or down, depending on the action of the controller. The integral action eliminates the steady error. The control parameters used in this work are presented in Table 9.

4. Results and Discussion

To demonstrate the model applicability, first, the model was validated under steady-state and transient behavior using manufacturer’s predictions [20], plant data [21, 22], and RELAP5 code [23]. Afterwards, the model was applied to the SBWR to analyze steady-state conditions at 100% of the rated power and to analyze the transient behavior during the closure of all main steam line isolation valves (MSIVs). Table 2 shows the nominal values of the SBWR used in the simulation.

4.1. Validation

The qualification assessment of the simplified model is presented in Appendix. The model predictions were compared with plant data and manufacturer’s predictions under steady-state behavior [20, 22]. To demonstrate the model applicability under transient behavior, an event occurred in the Laguna Verde nuclear power plant (LVNPP) was selected [21]. The plant data is compared with that obtained with RELAP5 [23] and the simplified model described in this work. The comparison results between RELAP5 and simplified model against plant data show that both computer codes have the capability to adequately predict the reactor behavior.

4.2. Steady-State Behavior

Table 2 shows the nominal values of the SBWR used in the simulation. The axial power distribution in the core (Figure 3) shows a high peak value in nodes 4 and 5 that correspond to the middle part of the fuel assembly. The lower power factor corresponds to the top and bottom of the assembly, where the neutron leakage is higher. This is a quite standard profile in a BWR, but using different schemes of the control rod the axial power shape is flatted along the fuel assembly in order to decrease the hot point in case of operational transient.

Figure 5 shows the behavior of the mixture density in the core. This is one of the most important parameters for the reactor that works under natural circulation, where the differences of the density of the coolant between the top and bottom of the core establish the natural flow.

The void fraction distribution (Figure 6) along the fuel assembly corresponds to a quiet standard BWR. Here, the heat transfer to the coolant increases the void fraction 70% before the dry step.

The radial and axial temperature distribution in the fuel is shown in Figure 7. The first radial node corresponds to the surface of the fuel cladding that is keeping about 560 K by the coolant flow. This temperature is the typical value during normal operation in the standard BWR. The node 8 corresponds to the center of the fuel pellet. The higher temperature corresponds to the nodes 4 and 5 because these nodes have the higher axial power factor (see Figure 3). Figure 7 shows that the natural circulation works efficiently to cool the surface temperature, far from the thermal limit even during transient conditions as show in Section 4.3.

The change in the enthalpy from subcooled liquid to saturated liquid during the natural circulation of the liquid phase is shown in Figure 8. This effect corresponds to the heating of the coolant during its travel along the core.

The transfer of energy from the fuel to the coolant due to heat transfer increases the enthalpy of the mixture node by node as shown in Figure 9. When 𝑓<𝑚, two-phase saturated exist, this takes place approximately at core node 3. If 𝑚<𝑓 and 𝑙<𝑓, liquid phase occurs in the core node 1, and subcooled boiling appears in the core nodes 2 and 3, respectively.

Figure 10 shows the convective heat transfer coefficient. The drastic increment of this parameter in the bottom of the core shows the importance of this mechanism to remove heat from the nuclear fuel during natural circulation. Specifically, in node 1, there is heat transfer in single liquid phase, and in nodes 2, 3, and 4 the heat transfer is governed by subcooled boiling. For the rest of the nodes the heat transfer regime is nucleate boiling.

The superficial velocity of the gas and liquid is shown in Figures 11 and 12, respectively. The superficial velocity of the gas increases due to the heating of the liquid following the curve of the void fraction. Therefore, the superficial velocity of the liquid tends to go down.

4.3. Closure of MISVs

The analysis of this transient assumes normal functioning of plant instrumentation and controls plant protection and reactor protection systems (RPSs). The sequence of the events of this transient is presented in Table 10.

The main steam isolation valves start to close in the first second with a delay of 3 seconds before closing completely, and the position switches on the valves initiate a reactor scram when the valves are less than 90% open. According to the behavior shown in Figures 13 and 14, the closure of these valves inhibits steam flow to the feed water turbines terminating feed water flow. Mitigation of pressure increase is accomplished by initiation of the reactor scram via MSIV position switches and the protection system (Figure 15) and the opening of the safety relief valves (SRVs) to limit system pressure (Figure 16). The peak of thermal power reaches 120% of rated power after approximately 3 seconds as shown in Figure 15, which is due to the drastic reduction of the void fraction in the core (Figure 17). This peak of power is limited by the set point scram at 118% of neutronic power. The peak pressure will still remain considerably below the ASME code limit of 1375 psig (9.48 MPa) as shown in Figure 18. The loss of feed water flow produces a reduction in the water lever in the reactor vessel as shown in Figure 19; the initiation of the emergency core cooling systems (ECCSs) occurs at approximately 33 seconds after detection; and the low level is not shown here. The scram of the reactor causes a reduction of the core flow (Figure 20). Figure 21 shows the heat transfer from the fuel to the coolant, where the shape follows the scram curve.

Figure 22 shows the behavior of the core flow in respect to void fraction. After the closure of the MISVs, the void fraction decreases suddenly due to pressure raise in the reactor vessel, which produces unbalance between the gravity head and the total loop pressure drop; therefore, the core flow decreases. This decrease can be observed in the figure, where finally the core flow reaches 70% of the rated core flow with 20% of the rated void fraction. This behavior is very important to core cooling. Figure 23 shows the neutronic power behavior in respect to core flow. In general terms, this figure shows that after shutdown the core flow keeps maintaining high values for low power values, which must guarantee the core cooling due to residual heat. The expanded graphic (Figure 23) shows that at the maximum power level the core flow increases slightly as expected.

5. Conclusions

A simplified model for transient analysis on boiling water reactor (BWR) working with natural circulation is described in this paper. This model was based on lumped and distributed parameter approximations, which includes vessel dome and downcomer, neutron kinetics, fuel temperature, lower and upper plena reactor core, and pressure and level controls.

The steady state and the transient of closure of main steam valves were presented and analyzed. The numerical results obtained with the model show a behavior physically consistent both in steady-state and transient behavior. Appendix A is a comparison of the simplified model predictions with plant data, manufacturer’s predictions, and RELAP5 under steady-state and transient conditions, showing good predictions between plant data and codes calculations. Therefore, the simplified model described in this paper provides reliable results to predict the behavior of a BWR working with natural circulation under normal and abnormal conditions.

Our results in the transient of closure of main steam valves show that the cooling system due to natural circulation in the SBWR is around 70% of the rated core flow. This can guarantee the core cooling due to residual heat.

Natural convection circulation of coolant in the reactor cooling system is used to some extent in all the reactor concepts, and many of these uses are innovative. Therefore, the need is to develop or confirm the design basis and to develop and qualify computer codes to enable reliable safety analysis.

Simplified model as shown in this work can be used to predict phenomena on SBWR behavior during the evaluation of the safety of a nuclear power plant which requires analysis of the plant’s responses to postulate equipment failures or malfunctions. Such analysis helps to determine the limit conditions for operation, limit safety system settings, and design specifications for safety-related components and systems to protect public health and safety.

Appendices

A. Qualification Assessment from Steady State

The profile of the void fraction prediction was compared with the manufacturer’s predictions [20]. The numerical results obtained with the simplified model show that void fraction behavior is similar to the results reported by the final safety analysis report [20] for a typical BWR. This assessment is shown in Figure 24. In this figure, it can be observed that manufacturer’s results underpredict the void fraction in the subcooled boiling region (approximated of 0.15 to 0.35 of core length). After this region, the profile is predicted with an accuracy of 0.323%.

428168.fig.024
Figure 24: Profile of the void fractions in comparison with manufacturer’s predictions [20].

The steam flow data design [22] in the reactor for a typical BWR is predicted with a maximum accuracy of 2.4%. This assessment is shown in Figure 25 for different load lines.

428168.fig.025
Figure 25: Steam flow comparison with data design [22].

B. Qualification Assessment of the Simplified Model under Transient Conditions

For the validation of the model under transient conditions, a transient occurred in the Laguna Verde nuclear power plant (LVNPP) was selected [21]. The plant data is compared with those obtained with RELAP5 [23] and the simplified model described in this work. This transient occurred in LVNPP Unit 2 on June 25th 1999, during the maintenance activities in the filters of the reactor water cleanup system (RWCU). A wrong sequence in the opening of valves of the RWCU caused a differential pressure and small pieces of the filter going up inside the reactor vessel. The signal of high radiation in the main steam lines causes a reactor SCRAM. This transient occurs at 90% of the full power.

Figure 26 shows the nodalization diagram of the LVNPP for the RELAP5 computer code. This is a typical boiling water reactor 5. In this model, the feed water, turbine, suppression pool, control rod drive system, and emergency cooling systems are represented by boundary conditions. Four fuel rod channels represent the core of the reactor [23].

428168.fig.026
Figure 26: Nodalization of Laguna Verde nuclear power plant with RELAP5 [21].

The signal of high radiation in the main steam lines causes the scram of the reactor and this is shown in Figure 27. This Figure shows a comparison among plant data, RELAP5, and simplified model results. The prediction of the simplified model is between RELAP5 calculation and plant data.

428168.fig.027
Figure 27: Thermal power behavior during the shutdown of the reactor.

The pressure of the reactor vessel is shown in Figure 28. The maximum peak of pressure according to the plant data is 75.69 Kg/cm2, and the simplified model shows an error of 2.5% below that value. RELAP5 calculation shows that the peak of pressure is well predicted. The opening of the safety relief valves causes a fast reduction of the pressure. In the simplified model, this reduction is slower than RELAP5 prediction.

428168.fig.028
Figure 28: Pressure rise of the reactor due to closure of main steam valves.

The core flow reduction due to the trip of the recirculation pumps by high pressure in the reactor vessel is shown in Figure 29. In this case, RELAP5 and the simplified model overpredicts the core flow, but during the first 5 seconds of the transient the results obtained with the simplified model are closer to the plant data, with a maximum difference of 12%.

428168.fig.029
Figure 29: Core flow reduction by recirculation pump trip.

Figure 30 shows the behavior of the steam flow. Again, the prediction of the simplified model is between plant data and RELAP5 calculation. The highest difference during the first 5 seconds of the transient is lower than 3%.

428168.fig.030
Figure 30: Steam flow reduction during the closure of main steam valves.

The comparison results between RELAP5 and simplified model against plant data shows that both computer codes have the capability to adequately predict the reactor behavior. One of the main advantages of the simplified model is that it spends less time to simulate different transients with a high reliability in results.

Nomenclature
Coefficients (see (24))
𝑎,𝑏,𝑐,𝑑:Fraction of wide open of the FW control valve
𝐴fw:Fraction of wide open of the MSL control valve
𝐴ms:Cross-sectional area 𝐴𝑥𝑠
(m2):Delayed neutron concentration
𝑐𝑖:Distribution parameter
𝐶𝑜:Single-phase friction factor
𝐶𝑓0:FW control valve-size coefficient 𝐶𝑣fw
(m2):MSL control valve-size coefficient 𝐶𝑣ms
(m2): Liquid specific heat 𝐶𝑝𝑙
(𝐽/kgK):Hydraulic diameter 𝐷
(m):Error signal
𝑒𝑝(𝑡):Convective boiling factor
𝐹:Axial power factor
𝐹(𝑧):Friction factor
𝑓:Mass flux 𝐺
(kg/sm2):Condensation parameter 𝐻𝑜
(1/sK):Convective heat transfer coefficient 𝐻
(J/sm2K):Enthalpy
(J/kg): Difference in specific enthalpy between saturated vapor and liquid 𝑓𝑔
(J/kg):Liquid enthalpy 𝑙
(J/kg): Vapor enthalpy 𝑣
(J/kg):Superficial velocity 𝑗
(m/s):Derivative control gain 𝐾𝐷
(s):Integral control gain 𝐾𝐼
(s1):Integral gain 𝐾𝐼𝑙
(s1):Proportional control gain
𝐾𝑝:Proportional gain
𝐾𝑝𝑙:Polarization
𝐾po:Support core plate loss coefficient 𝐾psn
(m4):Separador loss coefficient 𝐾sep
(m4): Liquid thermal conductivity 𝑘𝑙
(J/smK):Length in the upper plenum 𝐿ps
(m):Length in the steam separators 𝐿sep
(m): Mass of liquid 𝑚𝑙
(kg):Mass of vapor 𝑚𝑣
(kg):Normalized neutron flux
𝑛(𝑡):Instrument level 𝑁inst
(m):Liquid level in the vessel𝑁𝑙
(m):Normal level𝑁𝑜
(m):100% of rated power𝑃0
(J/s):Heated perimeter𝑃𝐻
(m): Reactor pressure𝑝Rx
(Pa):Pressure 𝑝
(Pa):Heat flux 𝑞
(J/sm2):Volumetric heat generation by gamma radiation𝑞
(J/sm3):Output signal
𝑆cp:Average fuel temperature𝑇𝑓
(K): Moderator temperature𝑇𝑚
(K):Saturation temperature𝑇𝑠
(K):Average drift velocity𝜈𝑔𝑗
(m/s):100% of rated of the mass flow𝑊𝑜
(kg/s): Core mass flow rate𝑊𝑐
(kg/s):Gas mass flow rate 𝑊gin
(kg/s):Liquid mass flow rate 𝑊𝑗𝑙
(kg/s): Vapor mass flow rate 𝑊𝑗𝑣
(kg/s): Liquid mass flow rate 𝑊lin
(kg/s):100% of mass flow rate 𝑊nom
(kg/s):Mass flow rate through the steam separator 𝑊sep
(kg/s):Vapor quality
𝑥: Martinelli factor.

Greek Symbols
𝑋𝑡𝑡:Neutron delay fraction
𝛽: Portion of neutrons generated by the 𝛽𝑖 group
𝑖th:Volumetric vapor generation rate Γ
(kg/m3)Core pressure drop Δ𝑝𝑐
(Pa):Drop pressure due to gravity Δ𝑝𝑔
(Pa):Time step Δ𝑡
(s):Cell length in axial direction Δ𝑧
(s):Void fraction
𝜀𝑔:Neutron generation time Λ
(s):Density 𝜌
(kg/m3): Core density in each node 𝜌𝑐
(kg/m3):Doppler effect 𝜌𝐷
(pcm):Downcomer density 𝜌dw
(kg/m3):Reactor control rods 𝜌cr
(kg/m3): Density of saturated liquid 𝜌𝑓
(kg/m3):Density of the saturated vapor 𝜌𝑔
(kg/m3):Rmoderador temperature reactivity 𝜌𝑚
(pcm):Steam-separator density 𝜌sep
(kg/m3): Void fraction reactivity ?𝑣
(pcm) Inertial term (?𝑙/𝐴)
(m1):Time constant of the FW enthalpy 𝑡fw𝑠
(s):Time constant of the MS valve 𝑡ms
(s):Time constant of the FW valve 𝑡fw
(s):Two-phase multiplier.

Subscripts
?2𝑓0,𝑗:Fuel, saturated liquid
𝑓:Gas phase
𝑔:Liquid phase
Mixture in two-phase flow.

Acknowledgments

The authors are grateful to Erick Gilberto Espinosa-Martínez for assistance provided during this study. Additional thanks are also due to the reviewers of the paper whose comments and observations helped to improve it substantially.