# SBWR Model for Steady-State

Authors: Gilberto Espinosa-Paredes; Alejandro Nuñez-Carrera
Journal: Science and Technology of Nuclear Installations
 (2008)

---

## 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.


---

## Body

## 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 Section2; 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 Figure1. 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.Figure 1
SBWR configuration and flow paths.Table1 shows the nominal values used in the simulation, which correspond to design
characteristics of the
SBWR [9].Table 1
Nominal parameters for SBWR [9].
ParameterUnitsThermal
power1800 MWthTotal
inlet core flow6666 kg/secPressure
at the dome of vessel7.07 MPaTotal
steam flow1100 kg/secFeed water
temperature488KActive
fuel length2.74mPressure
vessel internal height22m


## 3. Mathematical Model

Figure2 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.Table 2
SBWR submodels and approximations.
Mode*lApproximationVessel dome and
downcomerLumped parametersNatural circulation
loopLumped parametersNeutron
kineticsLumped parametersFuel rodDistributed parametersLower plenum,
upper plenum, and coreDistributed parameters*The whole system
includes pressure
control and level control.Figure 2
Schematic diagram of the SBWR.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 subscriptl) and gas- (denoted by subscript g) phase mass and energy balances, and a mixture momentum balance (denoted by subscript m): (1)∂∂t(ρgεg)+∂∂z(jgρg)=Γ,∂∂t(ρlεl)+∂∂z(jlρl)=−Γ,(2)∂∂t(ρlhlεl)−εl∂p∂t+∂∂z(ρlhljl)=ql′′PHAx−s+ql‴εl−Γhf,(3)∂∂t(ρmhm)−∂p∂t+Gm∂hm∂z=q′′PHAx−s+q‴,(4)∂Gm∂t+∂∂z(Gm2ρm)=−∂p∂z−fGm22ρmDh−gρm, where (5)ρm=ρgεg+ρlεl,(6)ρmhm=ρghgεg+ρlhlεl,(7)Gm=ρgjg+ρljl,(8)εg+εl=1.In these equations,ρ is density, εg is the void fraction; j is superficial velocity; Γ is the volumetric vapor generation rate; h is the enthalpy; p is the pressure; q″ is the heat flux; PH is the heated perimeter; Ax−s is the cross-sectional area; f is the friction factor; Dh is the hydraulic diameter; and q‴ is the volumetric heat generation by gamma
radiation. For unheated sections (e.g., downcomer, lower and upper plenums) q″=0.The equations of state are usually expressed
in terms of two independent variables in the formρl=ρl(p,hl) for density and Tl=Tl(p,hl) 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
ρk=ρk(p) to density and Tk=Tk(p) where k=l or g.
For the gas phase and in all cases hg=hg(p).Depending on the mixture enthalpyhm,
the following two cases were considered.Case 1.
Ifhf(p)<hm<hg(p),
the gas and liquid phases are saturated at reactor pressure and hl=hf(p).Case 2.
Ifhm<hf and  hl<hf,
liquid phase occurs first and,  subsequently, subcooled boiling appears.Here,hf and hg 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 ifhl=hf(p); ρl=ρf(p) in  (3) and combining it with (1), the
volumetric vapor generation rate Γ is given by: (9)Γ=1hfg[q′′PHAx−s+q‴+(1−ρgεg∂hg∂p−ρfεl∂hf∂p)∂p∂t−ρfjl∂hf∂z].In these equations,∂hg/∂z=0,
due that the gas phase is saturated at the vessel pressure, hfg is the difference in specific enthalpy of
saturated vapor and liquid. For subcooled boiling, the volumetric vapor
generation is given by [11]: (10)Γ=1hfg[q′′PHAx−s(1−hf−hlhf−hld)+q‴−Hohfgνfgεg(Ts−Tl)+εg(1−ρg∂hg∂p)∂p∂t],where Ho is the condensation parameter, and the liquid subcooling (hf−hld) at the initiation of subcooled boiling is
given by [12](11)hf−hld={0.0022q″DhCpl/kl,Pe<70,000,154q″/G,Pe>70,000.Here,kl is the liquid thermal conductivity; Cpl is the liquid specific heat; and Ts is the saturation temperature.The mass flow rate of each phase is,
respectively, given by(12)Wg=jgAx−sρg,Wl=jlAx−sρl,where jl=j−jg.
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)jg=(Coj+νgj)εg.The distribution parameterCo and the average drift velocity νgj 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: (14)j=[(1ρg−1ρl)Γ−(εgρg∂ρg∂p+εlρl∂ρl∂p)∂p∂t−ξsub]Δz+(Wginρg+Wlinρl)1Ax−s,where(15)ξsub={εlρl(∂ρl∂hl)∂hl∂t,Subcooledboiling,0,Thermodynamicequilibrium.Here, Δz is the cell length in axial direction; Wgin and Wlin 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 ∂p/∂t 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.Table 3
Thermal hydraulics parameters.
ParameterValueAx−s4.024m2*Dh0.0152mC01.2–0.2(ρg/ρl)0.5;    εg<0.65H00.24(s⋅oK)−1PH849.57mνgj1.53g0.2[σg(ρl−ρg)/ρl2]0.2;   εg<0.65*516 assemblies.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Δz.


### 3.2. Reactor Power Model

The reactor power is given by(16)P(t,z)=n(t)F(z)P0,where F(z) is the axial power factor (which is given by
Figure 3); P0 is nominal power; and n(t) is the normalized neutron flux, which is
calculated from a point reactor kinetics model with six groups of delayed
neutrons: (17)dn(t)dt=ρ(t)−βΛn(t)+∑i=16λici(t),(18)dci(t)dt=βiΛn(t)+λici(t),i=1,2,…,6,where ci is a delayed neutron concentration of the ith precursor group normalized with the
steady-state neutron density; ρ is the net reactivity; β is the neutron delay fraction; Λ is the neutron generation time; and βi is the portion of neutrons generated by the  ith  group. The initial conditions are
given by n(0)=n0 and ci(0)=βin0/Λλi.
The parameters of the kinetics model are presented in Table 4.Table 4
Neutron delay fraction and decay constant [13].
βiλiGroup 12.470×10−40.0127s−1Group 21.355×10−30.0317s−1Group 31.222×10−30.1150s−1Group 42.646×10−30.3110s−1Group 58.320×10−41.4000s−1Group 61.690×10−43.8700s−1β6.5×10−3Λ4.0×10−5s−1Figure 3
Axial power distribution.The net reactivity of the nuclear reactor includes
four main components: feedback reactivity due to the void fraction in the two-phase
flowρv,
Doppler effect ρD due to fuel temperature, moderator temperature ρm,
and reactor control rods ρcr.
Therefore, we express the feedback reactivities as: (19)ρ(t)=ρv(t)+ρD(t)+ρm(t)+ρcr(t).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.Table 5
Reactivity correlations.
Reactivity (pcm)Void
fractions reactivityρv(〈εg〉)=−4.0925−41.036〈εg〉+30.8223〈εg〉2.Doppler
effectρD=1.3−0.005〈Tf〉+2.0×10−6〈Tf〉2+−3.0×10−11〈Tf〉3−1.1×10−3〈εg〉(〈Tf〉−325.0).Moderator temperature reactivityρm=1.656−1.035×10−2〈Tm〉+4.0×10−5〈Tm〉2−7.33×10−8〈Tm3〉.Control rod drive
reactivityρcr=−11.235+57.73(1−CRD)1+0.937(1−CRD).


### 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 Figure4.Figure 4
Arrangement
of the computational cells of gap, clad, and fuel.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:(20)(ρCp)∂T∂t=k1r∂∂r(r∂T∂r)+q‴(t),atr0≤r≤r6,(21)I.C.T(r,0)=T(r),att=0,(22)B.C.1−k∂T∂r=H∞(T−Tm),atr=r6,(23)B.C.2∂T∂r=0,atr=r0.In these equations,r is the cylindrical radial coordinate; r0 and r6 are defined in Figure 4, q‴(t)=P(t)/Vf at each axial node, where P is given by (16); Tm is the moderator temperature; and H∞ 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:(24)ajTj−1t+Δt+bjTjt+Δt+cjTi,j+1t+Δt=dj, where Tj−1t+Δt,  Tjt+Δt,
and Tj+1t+Δt are unknowns; aj, bj, cj,
and dj are coefficients, which are computed at the time t.
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 aj, bj, cj
are  dependents of thermophysical properties, that is, thermal conductivity,
density, and specific heat, and since they are at function of Tjt+Δt,
at least one iteration is needed.The interaction of the model with the
neutron kinetics model is through evaluation of the average fuel temperature〈Tf〉 and the volumetric heat rate generation q‴.
The interaction with the nuclear thermalhydraulics is through evaluation of the
coolant temperature and heat transfer coefficient. The average fuel temperature
is given by: (25)〈Tf〉=1Vf∫VfTfdV≃1Vf∑i=1n(V2Ti,2+V3Ti,3+V4Ti,4+V5Ti,5), where Vf=V2+V3+V4+V5  is the fuel volume. Here 2,3,4,  and 5
represent  the radial node according to Figure 4, while that i represents the axial index of the node. Fuel
parameters are shown in Figure 4 and Table 6.Table 6
Fuel model parameters [14].
FuelKf=3825.02〈Tf〉+129.41+6.08×10−11〈Tf〉3,J/m⋅s,(ρCp)f=8.5103×1011e535.28/〈Tf〉〈Tf〉(e535.28/〈Tf〉−1)2+2.4348×102〈Tf〉+1.6609×1016e−1.8970x104/〈Tf〉〈Tf〉2,J/m3⋅K,Vf=8.9335m3.CladKcl=7.51−0.0209Tcl−1.45×10−5Tcl2+7.67×10−9Tcl3,J/m⋅s,(ρCp)cl=1.8204×106+3.0386×105θcl−1.0637×105θcl2+2.8102×104θcl3+−2723.6θcl4,J/m3⋅K,θcl=Tcl−300200.GapKga=Hgapδgap=2936.81.14×10−4,J/m⋅s.


### 3.4. Convective Heat Transfer
Coefficients

The convective heat transfer coefficient for the single-phase liquid
and vapor (H1ϕ) is obtained from the classical
Dittus-Boelter correlation [16]. The heat transfer coefficient in both
subcooled boiling and nucleate boiling regimes (H2ϕ)  is calculated with the Chen correlation [17].
The convective heat transfer coefficient (H∞) is given by the maximum between H1ϕ and H2ϕ to avoid discontinuities, that is, H∞=MAX(H1ϕ,H2ϕ).The interaction of the convective heat transfer coefficients with the fuel heat transfer model is through the evaluation of the
clad temperatureT∣r=r6 (Figure 4). The convective heat transfer
coefficient H∞ 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)dpdt=−νl∑jWjl+(∂νl/∂hl)p∑jWjl(hjl−hl)Q+Q′+νv∑jWjv+(∂νv/∂hv)p∑jWjv(hjv−hv)Q+Q′, where Q denotes ml[νl(∂νl/∂hl)p+(∂νl/∂p)hl]   and Q′ denotes mv[νv(∂νv/∂hv)p+(∂νv/∂p)hv], and    Wjl and Wjv are the mass flow rates into or out of the
liquid and vapor regions, respectively, and
hjl and hjv are the enthalpy of the fluid entering or
leaving the liquid and vapor regions, respectively. The mass of liquid ml,
mass of vapor mv,
enthalpy of liquid hl,
and enthalpy of vapor hv are given by the following balance equations: (27)dmldt=∑jlWjl,dmvdt=∑jvWjv,(28)dhldt=1ml∑jlWjl(hjl−hl)+νldpdt,dhvdt=1mv∑jvWjv(hjv−hv)+νvdpdt.The flashing and condensation flows are
given by: (29)Wf=mlΔt(hl−hfhg−hf),Wco=mvΔt(hg−hvhg−hf).Here Δt is the time step. The density and temperature
in the downcomer is given by: (30)ρdw=ρdw(p,hl),(31)Tdw=Tdw(p,hl).The level in the vessel is calculated as a
function of liquid volume, that is,Nl=Nl(Vl),
where the liquid volume is given by Vl=mlνl.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
Figure1. 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: (32)dWcdt=(∑i=1nliAi)−1(−KpsnWc2ρdw−KsepWsep2ρsep−Δpc+Δpg), where (∑l/A) is the inertial term; ρdw is the downcomer density (30); ρsep is the steam-separator density; Wsep is the flow mass through the steam separator; Kpsn is the support core plate loss coefficient; Ksep is the separator loss coefficient; Δpc is the core pressure drop; and Δpg is the pressure drop due to gravity.The total
core pressure drop is the sum of the frictional, acceleration, and
gravitational components:(33)−Δpc=∑j=112[ϕf0,j22Cf0Wc2ΔzAx−s2Dhρf+Wc2Ax−s2(1ρc,j+1−1ρc,j)+ρc,jgΔz],where  j=1,2,…,12 are used to indicate the
core node number; Cf0 is the single-phase friction factor; Δzj is the core node length; Dh is the hydraulic diameter; ρf is the density of saturated liquid; ρc is the core density in each node, which is
given by (5); and ϕf0,j2 is the two-phase multiplier given by
[18]: (34)ϕf02=1+x(ρfρg−1),where x  is the vapor quality and ρg is the density of the saturated vapor.The pressure drop due to gravity
without core effects, which are considered in (32), is given by: (35)Δpg=gNl(Vl)ρdw−g(Lps+Lsep)ρsep, where Nl, Lps,
and Lsep 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.Table 7
Circulation natural parameters.
ParameterValue(∑l/A)9m−1Kpsn1.3m−4Ksep0.01m−4Lps8mLsep5mCf00.005Δz0.2283mThe 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]: (36)Wms=Cvmsf(Ams*)Δpmsρv|Δpms|,where Cvms=0.4158 is the control valve size coefficient, whose
value includes the MSL system; f(Ams*) is the fraction of the total flow area; Ams* is the fraction of wide open; Δpms=pRx−6.7203×106 is the MSL pressure drop; pRx 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:(37)Wfw=Cvfwg(Afw*)ρlΔpfw,where Cvfw=0.04595 is the control valve size coefficient, whose
value includes the FW system; Δpfw=7.9976×106Pa−pRx 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]: (38)haa=−187903.0+50281.0(Wms*)−983.3(Wms*)2+8.88(Wms*)3−0.03(Wms*)4,where Wms*=(Wms/Wnom)×100, Wms is given by (43), and Wnom is 100% of mass flow rate. The dynamic effects
are considered as a response of first-order system [19]: (39)dhfwdt=haa−hfwτfws, where τfws is the FW time constant. Table 8 presents the
time constant.Table 8
Time constant used in this work.
τfws1000.0sτms1.25sτfw10.0sThe dynamic of the valves was modeled as a
response of a first-order system:(40)dAms*dt=Scp−Ams*τms,dAfw*dt=Scl−Afw*τfw,where τms and τfw are time constants of the MS valve and FW
valve, respectively; Scp
and
Scl 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 signalScp in direct proportion to the error signal ep(t) [19]: (41)Scp=Kpo+Kpep(t)+KDdeldt+KI∫ep(t)dt, where Kpo is the polarization and the controller value
output when there is no error; Kp, KD,
and KI are, respectively, the proportional,
derivative, and integral control gains. The tracking error is given by: (42)ep(t)=pRx−7.07×106Papnom. Here, pRx is the reactor pressure; 7.07×106Pa   is the set point of pressure; and pnom is the nominal pressure. The proportional-derivative feedback control action Kpep+KD(del/dt) is used to induce certain stability margin in the controlled pressure. However, the stability of pressure does not ensure that the tracking error ep(t) can be made as small as desired. To this end,
the integral action KI∫el(t)dt,
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)Scl(t)=Kplel(t)+KIl∫el(t)dt, where Kpl is the proportional gain; KIl is the integral gain; and the tracking error is given by: (44)el(t)=Wms−WfwWo+Ninst−0.97mNo, where Wo is the mass flow rate; No is the normal level; Ninst 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.Table 9
Control parameters.
Kpo0.75Kp18.4KI18.4Kpl0.90KIl0.02


## 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 subscriptl) and gas- (denoted by subscript g) phase mass and energy balances, and a mixture momentum balance (denoted by subscript m): (1)∂∂t(ρgεg)+∂∂z(jgρg)=Γ,∂∂t(ρlεl)+∂∂z(jlρl)=−Γ,(2)∂∂t(ρlhlεl)−εl∂p∂t+∂∂z(ρlhljl)=ql′′PHAx−s+ql‴εl−Γhf,(3)∂∂t(ρmhm)−∂p∂t+Gm∂hm∂z=q′′PHAx−s+q‴,(4)∂Gm∂t+∂∂z(Gm2ρm)=−∂p∂z−fGm22ρmDh−gρm, where (5)ρm=ρgεg+ρlεl,(6)ρmhm=ρghgεg+ρlhlεl,(7)Gm=ρgjg+ρljl,(8)εg+εl=1.In these equations,ρ is density, εg is the void fraction; j is superficial velocity; Γ is the volumetric vapor generation rate; h is the enthalpy; p is the pressure; q″ is the heat flux; PH is the heated perimeter; Ax−s is the cross-sectional area; f is the friction factor; Dh is the hydraulic diameter; and q‴ is the volumetric heat generation by gamma
radiation. For unheated sections (e.g., downcomer, lower and upper plenums) q″=0.The equations of state are usually expressed
in terms of two independent variables in the formρl=ρl(p,hl) for density and Tl=Tl(p,hl) 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
ρk=ρk(p) to density and Tk=Tk(p) where k=l or g.
For the gas phase and in all cases hg=hg(p).Depending on the mixture enthalpyhm,
the following two cases were considered.Case 1.
Ifhf(p)<hm<hg(p),
the gas and liquid phases are saturated at reactor pressure and hl=hf(p).Case 2.
Ifhm<hf and  hl<hf,
liquid phase occurs first and,  subsequently, subcooled boiling appears.Here,hf and hg 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 ifhl=hf(p); ρl=ρf(p) in  (3) and combining it with (1), the
volumetric vapor generation rate Γ is given by: (9)Γ=1hfg[q′′PHAx−s+q‴+(1−ρgεg∂hg∂p−ρfεl∂hf∂p)∂p∂t−ρfjl∂hf∂z].In these equations,∂hg/∂z=0,
due that the gas phase is saturated at the vessel pressure, hfg is the difference in specific enthalpy of
saturated vapor and liquid. For subcooled boiling, the volumetric vapor
generation is given by [11]: (10)Γ=1hfg[q′′PHAx−s(1−hf−hlhf−hld)+q‴−Hohfgνfgεg(Ts−Tl)+εg(1−ρg∂hg∂p)∂p∂t],where Ho is the condensation parameter, and the liquid subcooling (hf−hld) at the initiation of subcooled boiling is
given by [12](11)hf−hld={0.0022q″DhCpl/kl,Pe<70,000,154q″/G,Pe>70,000.Here,kl is the liquid thermal conductivity; Cpl is the liquid specific heat; and Ts is the saturation temperature.The mass flow rate of each phase is,
respectively, given by(12)Wg=jgAx−sρg,Wl=jlAx−sρl,where jl=j−jg.
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)jg=(Coj+νgj)εg.The distribution parameterCo and the average drift velocity νgj 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: (14)j=[(1ρg−1ρl)Γ−(εgρg∂ρg∂p+εlρl∂ρl∂p)∂p∂t−ξsub]Δz+(Wginρg+Wlinρl)1Ax−s,where(15)ξsub={εlρl(∂ρl∂hl)∂hl∂t,Subcooledboiling,0,Thermodynamicequilibrium.Here, Δz is the cell length in axial direction; Wgin and Wlin 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 ∂p/∂t 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.Table 3
Thermal hydraulics parameters.
ParameterValueAx−s4.024m2*Dh0.0152mC01.2–0.2(ρg/ρl)0.5;    εg<0.65H00.24(s⋅oK)−1PH849.57mνgj1.53g0.2[σg(ρl−ρg)/ρl2]0.2;   εg<0.65*516 assemblies.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Δz.


## 3.2. Reactor Power Model

The reactor power is given by(16)P(t,z)=n(t)F(z)P0,where F(z) is the axial power factor (which is given by
Figure 3); P0 is nominal power; and n(t) is the normalized neutron flux, which is
calculated from a point reactor kinetics model with six groups of delayed
neutrons: (17)dn(t)dt=ρ(t)−βΛn(t)+∑i=16λici(t),(18)dci(t)dt=βiΛn(t)+λici(t),i=1,2,…,6,where ci is a delayed neutron concentration of the ith precursor group normalized with the
steady-state neutron density; ρ is the net reactivity; β is the neutron delay fraction; Λ is the neutron generation time; and βi is the portion of neutrons generated by the  ith  group. The initial conditions are
given by n(0)=n0 and ci(0)=βin0/Λλi.
The parameters of the kinetics model are presented in Table 4.Table 4
Neutron delay fraction and decay constant [13].
βiλiGroup 12.470×10−40.0127s−1Group 21.355×10−30.0317s−1Group 31.222×10−30.1150s−1Group 42.646×10−30.3110s−1Group 58.320×10−41.4000s−1Group 61.690×10−43.8700s−1β6.5×10−3Λ4.0×10−5s−1Figure 3
Axial power distribution.The net reactivity of the nuclear reactor includes
four main components: feedback reactivity due to the void fraction in the two-phase
flowρv,
Doppler effect ρD due to fuel temperature, moderator temperature ρm,
and reactor control rods ρcr.
Therefore, we express the feedback reactivities as: (19)ρ(t)=ρv(t)+ρD(t)+ρm(t)+ρcr(t).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.Table 5
Reactivity correlations.
Reactivity (pcm)Void
fractions reactivityρv(〈εg〉)=−4.0925−41.036〈εg〉+30.8223〈εg〉2.Doppler
effectρD=1.3−0.005〈Tf〉+2.0×10−6〈Tf〉2+−3.0×10−11〈Tf〉3−1.1×10−3〈εg〉(〈Tf〉−325.0).Moderator temperature reactivityρm=1.656−1.035×10−2〈Tm〉+4.0×10−5〈Tm〉2−7.33×10−8〈Tm3〉.Control rod drive
reactivityρcr=−11.235+57.73(1−CRD)1+0.937(1−CRD).


## 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 Figure4.Figure 4
Arrangement
of the computational cells of gap, clad, and fuel.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:(20)(ρCp)∂T∂t=k1r∂∂r(r∂T∂r)+q‴(t),atr0≤r≤r6,(21)I.C.T(r,0)=T(r),att=0,(22)B.C.1−k∂T∂r=H∞(T−Tm),atr=r6,(23)B.C.2∂T∂r=0,atr=r0.In these equations,r is the cylindrical radial coordinate; r0 and r6 are defined in Figure 4, q‴(t)=P(t)/Vf at each axial node, where P is given by (16); Tm is the moderator temperature; and H∞ 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:(24)ajTj−1t+Δt+bjTjt+Δt+cjTi,j+1t+Δt=dj, where Tj−1t+Δt,  Tjt+Δt,
and Tj+1t+Δt are unknowns; aj, bj, cj,
and dj are coefficients, which are computed at the time t.
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 aj, bj, cj
are  dependents of thermophysical properties, that is, thermal conductivity,
density, and specific heat, and since they are at function of Tjt+Δt,
at least one iteration is needed.The interaction of the model with the
neutron kinetics model is through evaluation of the average fuel temperature〈Tf〉 and the volumetric heat rate generation q‴.
The interaction with the nuclear thermalhydraulics is through evaluation of the
coolant temperature and heat transfer coefficient. The average fuel temperature
is given by: (25)〈Tf〉=1Vf∫VfTfdV≃1Vf∑i=1n(V2Ti,2+V3Ti,3+V4Ti,4+V5Ti,5), where Vf=V2+V3+V4+V5  is the fuel volume. Here 2,3,4,  and 5
represent  the radial node according to Figure 4, while that i represents the axial index of the node. Fuel
parameters are shown in Figure 4 and Table 6.Table 6
Fuel model parameters [14].
FuelKf=3825.02〈Tf〉+129.41+6.08×10−11〈Tf〉3,J/m⋅s,(ρCp)f=8.5103×1011e535.28/〈Tf〉〈Tf〉(e535.28/〈Tf〉−1)2+2.4348×102〈Tf〉+1.6609×1016e−1.8970x104/〈Tf〉〈Tf〉2,J/m3⋅K,Vf=8.9335m3.CladKcl=7.51−0.0209Tcl−1.45×10−5Tcl2+7.67×10−9Tcl3,J/m⋅s,(ρCp)cl=1.8204×106+3.0386×105θcl−1.0637×105θcl2+2.8102×104θcl3+−2723.6θcl4,J/m3⋅K,θcl=Tcl−300200.GapKga=Hgapδgap=2936.81.14×10−4,J/m⋅s.


## 3.4. Convective Heat Transfer
Coefficients

The convective heat transfer coefficient for the single-phase liquid
and vapor (H1ϕ) is obtained from the classical
Dittus-Boelter correlation [16]. The heat transfer coefficient in both
subcooled boiling and nucleate boiling regimes (H2ϕ)  is calculated with the Chen correlation [17].
The convective heat transfer coefficient (H∞) is given by the maximum between H1ϕ and H2ϕ to avoid discontinuities, that is, H∞=MAX(H1ϕ,H2ϕ).The interaction of the convective heat transfer coefficients with the fuel heat transfer model is through the evaluation of the
clad temperatureT∣r=r6 (Figure 4). The convective heat transfer
coefficient H∞ 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)dpdt=−νl∑jWjl+(∂νl/∂hl)p∑jWjl(hjl−hl)Q+Q′+νv∑jWjv+(∂νv/∂hv)p∑jWjv(hjv−hv)Q+Q′, where Q denotes ml[νl(∂νl/∂hl)p+(∂νl/∂p)hl]   and Q′ denotes mv[νv(∂νv/∂hv)p+(∂νv/∂p)hv], and    Wjl and Wjv are the mass flow rates into or out of the
liquid and vapor regions, respectively, and
hjl and hjv are the enthalpy of the fluid entering or
leaving the liquid and vapor regions, respectively. The mass of liquid ml,
mass of vapor mv,
enthalpy of liquid hl,
and enthalpy of vapor hv are given by the following balance equations: (27)dmldt=∑jlWjl,dmvdt=∑jvWjv,(28)dhldt=1ml∑jlWjl(hjl−hl)+νldpdt,dhvdt=1mv∑jvWjv(hjv−hv)+νvdpdt.The flashing and condensation flows are
given by: (29)Wf=mlΔt(hl−hfhg−hf),Wco=mvΔt(hg−hvhg−hf).Here Δt is the time step. The density and temperature
in the downcomer is given by: (30)ρdw=ρdw(p,hl),(31)Tdw=Tdw(p,hl).The level in the vessel is calculated as a
function of liquid volume, that is,Nl=Nl(Vl),
where the liquid volume is given by Vl=mlνl.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
Figure1. 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: (32)dWcdt=(∑i=1nliAi)−1(−KpsnWc2ρdw−KsepWsep2ρsep−Δpc+Δpg), where (∑l/A) is the inertial term; ρdw is the downcomer density (30); ρsep is the steam-separator density; Wsep is the flow mass through the steam separator; Kpsn is the support core plate loss coefficient; Ksep is the separator loss coefficient; Δpc is the core pressure drop; and Δpg is the pressure drop due to gravity.The total
core pressure drop is the sum of the frictional, acceleration, and
gravitational components:(33)−Δpc=∑j=112[ϕf0,j22Cf0Wc2ΔzAx−s2Dhρf+Wc2Ax−s2(1ρc,j+1−1ρc,j)+ρc,jgΔz],where  j=1,2,…,12 are used to indicate the
core node number; Cf0 is the single-phase friction factor; Δzj is the core node length; Dh is the hydraulic diameter; ρf is the density of saturated liquid; ρc is the core density in each node, which is
given by (5); and ϕf0,j2 is the two-phase multiplier given by
[18]: (34)ϕf02=1+x(ρfρg−1),where x  is the vapor quality and ρg is the density of the saturated vapor.The pressure drop due to gravity
without core effects, which are considered in (32), is given by: (35)Δpg=gNl(Vl)ρdw−g(Lps+Lsep)ρsep, where Nl, Lps,
and Lsep 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.Table 7
Circulation natural parameters.
ParameterValue(∑l/A)9m−1Kpsn1.3m−4Ksep0.01m−4Lps8mLsep5mCf00.005Δz0.2283mThe 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]: (36)Wms=Cvmsf(Ams*)Δpmsρv|Δpms|,where Cvms=0.4158 is the control valve size coefficient, whose
value includes the MSL system; f(Ams*) is the fraction of the total flow area; Ams* is the fraction of wide open; Δpms=pRx−6.7203×106 is the MSL pressure drop; pRx 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:(37)Wfw=Cvfwg(Afw*)ρlΔpfw,where Cvfw=0.04595 is the control valve size coefficient, whose
value includes the FW system; Δpfw=7.9976×106Pa−pRx 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]: (38)haa=−187903.0+50281.0(Wms*)−983.3(Wms*)2+8.88(Wms*)3−0.03(Wms*)4,where Wms*=(Wms/Wnom)×100, Wms is given by (43), and Wnom is 100% of mass flow rate. The dynamic effects
are considered as a response of first-order system [19]: (39)dhfwdt=haa−hfwτfws, where τfws is the FW time constant. Table 8 presents the
time constant.Table 8
Time constant used in this work.
τfws1000.0sτms1.25sτfw10.0sThe dynamic of the valves was modeled as a
response of a first-order system:(40)dAms*dt=Scp−Ams*τms,dAfw*dt=Scl−Afw*τfw,where τms and τfw are time constants of the MS valve and FW
valve, respectively; Scp
and
Scl 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 signalScp in direct proportion to the error signal ep(t) [19]: (41)Scp=Kpo+Kpep(t)+KDdeldt+KI∫ep(t)dt, where Kpo is the polarization and the controller value
output when there is no error; Kp, KD,
and KI are, respectively, the proportional,
derivative, and integral control gains. The tracking error is given by: (42)ep(t)=pRx−7.07×106Papnom. Here, pRx is the reactor pressure; 7.07×106Pa   is the set point of pressure; and pnom is the nominal pressure. The proportional-derivative feedback control action Kpep+KD(del/dt) is used to induce certain stability margin in the controlled pressure. However, the stability of pressure does not ensure that the tracking error ep(t) can be made as small as desired. To this end,
the integral action KI∫el(t)dt,
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)Scl(t)=Kplel(t)+KIl∫el(t)dt, where Kpl is the proportional gain; KIl is the integral gain; and the tracking error is given by: (44)el(t)=Wms−WfwWo+Ninst−0.97mNo, where Wo is the mass flow rate; No is the normal level; Ninst 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.Table 9
Control parameters.
Kpo0.75Kp18.4KI18.4Kpl0.90KIl0.02


## 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

Table2 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.Figure5 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.Figure 5
Axial distribution of the mixture
density.The void fraction
distribution (Figure6) 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.Figure 6
Axial distribution of the void
fraction.The radial and axial temperature
distribution in the fuel is shown in Figure7. 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.Fuel rod temperature distribution. (a) Surface graphic, (b) contour lines.
(a)(b)The change in the
enthalpy from subcooled
liquid to saturated
liquid during the natural circulation of the liquid phase is shown in Figure8.
This effect corresponds to the heating of the coolant during its travel along
the core.Figure 8
Axial distribution of the enthalpy of
the liquid.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 Figure9. When hf<hm,
two-phase saturated exist, this takes place approximately at core node 3. If hm<hf and hl<hf,
liquid phase occurs in the core node 1, and subcooled boiling appears in the
core nodes 2 and 3, respectively.Figure 9
Axial distribution of the mixture.hf (= 1271026) liquid saturated enthalpy.Figure10 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.Figure 10
Axial distribution of the convective
heat transfer during natural circulation.The superficial velocity
of the gas and liquid is
shown in Figures11 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.Figure 11
Superficial velocity for the gas phase.Figure 12
Superficial velocity for the liquid
phase.


### 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 Table10.Table 10
Sequence of events of closure of main steam isolation valves.
Elapsed
time (sec.)Event0.0Reactor under steady-state
condition rate at full power.1.0Closure of all main
steam isolation valves.2.0Maximum power peak is reached (120% of rated).2.0The peak of pressure
is maximum, and the SRVs opening to relief pressure.2.0+Reactor lever decreases until the emergency
core cooling system actuation.10.0+Core flow tend to stabilize at new equilibrium
conditions
(natural circulation).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 Figures13 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 13
Feed water behavior.Figure 14
Steam flow behavior.Figure 15
Neutronic power behavior during the
closure of MSIVs.Figure 16
SRV'S steam flow by closure of MSIVs.Figure 17
Average void fraction behavior.Figure 18
Pressure change in the dome of the
reactor vessel.Figure 19
Water lever reduction due to loss of
feed water.Figure 20
Core flow behavior.Figure 21
Heat flux behavior.Figure22 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.Figure 22
Core flow versus void fraction.Figure 23
Neutronic power versus core flow.


## 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

Table2 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.Figure5 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.Figure 5
Axial distribution of the mixture
density.The void fraction
distribution (Figure6) 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.Figure 6
Axial distribution of the void
fraction.The radial and axial temperature
distribution in the fuel is shown in Figure7. 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.Fuel rod temperature distribution. (a) Surface graphic, (b) contour lines.
(a)(b)The change in the
enthalpy from subcooled
liquid to saturated
liquid during the natural circulation of the liquid phase is shown in Figure8.
This effect corresponds to the heating of the coolant during its travel along
the core.Figure 8
Axial distribution of the enthalpy of
the liquid.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 Figure9. When hf<hm,
two-phase saturated exist, this takes place approximately at core node 3. If hm<hf and hl<hf,
liquid phase occurs in the core node 1, and subcooled boiling appears in the
core nodes 2 and 3, respectively.Figure 9
Axial distribution of the mixture.hf (= 1271026) liquid saturated enthalpy.Figure10 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.Figure 10
Axial distribution of the convective
heat transfer during natural circulation.The superficial velocity
of the gas and liquid is
shown in Figures11 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.Figure 11
Superficial velocity for the gas phase.Figure 12
Superficial velocity for the liquid
phase.


## 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 Table10.Table 10
Sequence of events of closure of main steam isolation valves.
Elapsed
time (sec.)Event0.0Reactor under steady-state
condition rate at full power.1.0Closure of all main
steam isolation valves.2.0Maximum power peak is reached (120% of rated).2.0The peak of pressure
is maximum, and the SRVs opening to relief pressure.2.0+Reactor lever decreases until the emergency
core cooling system actuation.10.0+Core flow tend to stabilize at new equilibrium
conditions
(natural circulation).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 Figures13 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 13
Feed water behavior.Figure 14
Steam flow behavior.Figure 15
Neutronic power behavior during the
closure of MSIVs.Figure 16
SRV'S steam flow by closure of MSIVs.Figure 17
Average void fraction behavior.Figure 18
Pressure change in the dome of the
reactor vessel.Figure 19
Water lever reduction due to loss of
feed water.Figure 20
Core flow behavior.Figure 21
Heat flux behavior.Figure22 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.Figure 22
Core flow versus void fraction.Figure 23
Neutronic power versus core flow.


## 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.



---

*Source: 428168-2008-05-11.xml*