Numerical Simulations and Design Optimization of the PHT Loop of Natural Circulation BWR

Mathematical modeling and numerical simulation of natural circulation boiling water reactor (NCBWR) are very important in order to study its performance for different designs and various off-design conditions and for design optimization. In the present work, parametric studies of the primary heat transport loop of NCBWR have been performed using lumped parameter models and RELAP5/MOD3.4 code. The lumped parameter models are based on the drift flux model and homogeneous equilibrium mixture (HEM) model of two-phase flow. Numerical simulations are performed with both models. Compared to the results obtained from the HEM model, those obtained from the drift flux model are closer to RELAP5. The variations of critical heat flux with various geometric parameters and operating conditions are thoroughly investigated. The material required to construct the primary heat transport (PHT) loop of NCBWR has been minimized using sequential quadratic programming. The stability of NCBWR has also been verified at the optimum point.


INTRODUCTION
Water-cooled nuclear power plants (light and heavy water reactors) represent 94% of the global nuclear power capacity, and advanced designs are being developed to meet future energy needs.Common goals for advanced designs are high reliability and competitive economics while meeting stringent safety requirements.Passive safety systems based on natural circulation are the key to several evolutionary water-cooled designs and many revolutionary water-cooled reactor designs.In natural circulation boiling water reactor (NCBWR), the heat removal from the core takes place by natural circulation during startup, power rising, and accidental conditions in addition to the rated full power operating conditions.The thermal hydraulic analysis of NCBWR is very important for better insight into the process that governs the in-vessel physics.Parametric studies of NCBWR will guide the designer to choose the best possible design.Numerical simulations and parametric studies can be performed using thermal hydraulic codes such as RELAP5.This code is based on governing equations of thermal hydraulics and takes considerable computational time.Extensive parametric studies and design optimization can be done with the help of computationally inexpensive mathematical model, such as lumped parameter model (LPM).
Aybar et al. [1] have investigated the performance of Ohio State University inherently safe reactor (OSU-ISR) as a function of several design variables, to know the interdependency among the several system variables.They simulated the steady-state performance of OSU-ISR with dynamic simulation for nuclear power plants (DSNP) code and validated the results with RELAP5.Jafari et al. [2] have analyzed the effect of design parameters on natural circulation loop and have done stability analysis using RELAP5 thermal hydraulic code.Rao et al. [3] have studied the performance of natural circulation loop using homogeneous equilibrium model (HEM) and drift flux model through one-dimensional approach.Kok and van der Hagen [4] developed a steady-state mathematical model to simulate DESIRE test facility and conducted parametric studies.Garg et al. [5] did parametric studies of NCBWR using artificial neural network (ANN) model, which was trained using the results of RELAP5/MOD3.2simulation.
They have illustrated a methodology of training ANN with RELAP, which could then be used for parametric studies and design optimization.
In the present study, steady-state lumped parameter models (LPMs) are derived for the primary heat transport (PHT) loop of NCBWR.The same loop is also modeled in RELAP5/MOD3.4code.Parametric trends are computed with both models (LPM and RELAP), and the results are compared.The variations of mass flow rate and void fraction with geometric properties and with the power are studied.After that the critical heat flux (CHF) phenomenon is investigated.The variation of CHF with various geometric and operating conditions is studied.Subsequently, design optimization of NCBWR is done to minimize the volume of material required to construct the PHT loop (referred to as the "minimum material condition" in this paper) subject to constraints based on safety considerations.The stability of the NCBWR is checked at minimum material condition, and alternative designs are suggested based on stability considerations.

Description of the system
In natural circulation,the primary pumps are eliminated, and the necessary flow rates are achieved by locating the steam drums at a suitable height from the center of the core, taking advantage of the reactor building height.A schematic view of a natural circulation boiling water reactor (NCBWR) is shown in Figure 1.The input data are listed in Tables 1 and  2. These data are for a typical design of NCBWR, taken from literature [6].The PHT loop consists of core, tail pipes (risers), steam drum, downcomer, inlet header, and inlet feed water pipes.Constant heat is supplied in the active core region.The steam water mixture from each coolant channel is led through tail pipes to the steam drum, which is located at certain elevation with respect to inlet header (coolant channel bottom).In the present work, 113 coolant channels are taken for the analysis.
The two-phase mixture is present in the core and riser sections, and the single phase liquid is in the downcomer section.Due to density difference between these two sections natural circulation will develop.The natural circulation flow rate is decided by the density difference of the fluid between the downcomer (cold leg) and that of the riser (hot leg), and the flow resistances.The enthalpy of the coolant at the core inlet depends on the feed water flow rate, its temperature, and the recirculation flow rate.

RELAP modeling
The PHT loop is modeled using RELAP5/MOD3.4.It is a time domain code, which is based on two-fluid formulation.The nodalization scheme for the PHT loop is shown in Figure 2. The PHT loop components are modeled in RELAP using pipe components.The boundary conditions, that is, steam outlet from steam drum and the feed water inlet are defined using time dependent volumes.The pipe component  P1 is divided into 12 volumes.The first two volumes model the feed water inlet pipe, the third volume models the bottom of the core followed by 5 volumes of active core, one volume of core top and three volumes for tail pipes.The steam drum (pipe component P2) is divided into 10 volumes.The saturated water from steam drum and the feed water are mixed in single volume S1.The downcomer is divided into 2 volumes, and the inlet header is modeled as a single volume.
All the components are connected using single junctions.The heat supply to the active core is modeled using a heat structure component divided into 5 nodes (equivalent to 5 volumes of active core), and constant power option (no neutronics) is used.

Steady-state modeling
Two lumped parametermodels are derived using drift flux and homogeneous equilibrium model (HEM) from the onedimensional steady-state equations of mass, momentum, and energy balances.The following assumptions are made while deriving the equations.These assumptions are used to reduce the complexity of the system at the expense of some accuracy and to give a better insight into the system.
(1) The thermodynamic properties over the entire loop are evaluated at the system pressure.
(2) The liquid and vapor phases are in thermodynamic equilibrium with each other.
(3) The processes of subcooled boiling and flashing are not considered.
(4) Viscous dissipation and changes in kinetic energy, potential energy, and flow work are neglected in the energy equation.
(5) The distribution of linear power of the fuel in the axial direction is uniform.
(6) The coolant in the downcomer is in the liquid phase only.
The present models are similar to an earlier lumped parameter dynamic model of natural circulation loop derived by van Bragt and van der Hagen [7].The core is divided into single-phase region (from core inlet to boiling boundary) and two-phase region (from boiling boundary to the core exit).The one-dimensional unsteady-state continuity, momentum, and energy equations are given below: (1) The steady-state equations are obtained by equating the time derivative terms to zero in the above three equations.
The boiling boundary equation is derived by integrating the energy equation from the core inlet to the boiling boundary: The expression for the exit steam quality is derived by integrating the energy equation from boiling boundary to exit of the core: Under steady-state conditions, the core exit mass flux density is equal to core inlet mass flux density.Applying steady-state mass and energy balance at the junction of steam drum, feedwater and downcomer, the core inlet enthalpy can be expressed in terms of feedwater inlet enthalpy as The relation between volumetric quality (β) and void fraction (α) is given by Zuber and Findlay [8] as follows: where C 0 is the flow distribution parameter, V g j is drift velocity, and j is the volumetric flux.The volumetric quality β can be expressed in terms of core exit quality as For the drift flux model, numerous empirical correlations are available for C 0 and V g j in literature.Coddington and Macian [9] discussed 13 correlations and their ranges of operation.They found that the void fraction predictions with Maier and Coddington correlation are better than those for the other correlations.Therefore, in the present study, Maier and Coddington correlation has been used.This correlation was purely based on empirical approach to the evaluation of C 0 and V g j of the drift-flux equation.The distribution parameter is a function of pressure, and drift velocity is a function of pressure and mass flux density.Under steadystate and constant pressure conditions and with constant linear power, the axial variation of mass quality in the core will be linear.Using the above relations for void fraction α and volumetric quality β, the expression for the core average void fraction is obtained as where Under steady-state and constant pressure conditions, the riser void fraction is equal to the exit void fraction of the core.The momentum equation is integrated over the natural circulation loop, including the core, the riser, and the downcomer.Since the loop is closed, the total pressure drop in the natural circulation loop must be zero.There are various pressure drops in different sections of the loop, arising due to gravitational, frictional, acceleration, and minor losses.These pressure drop terms are derived by integrating the respective terms of the one-dimensional momentum equation along each component of the boiling system.The friction factor is calculated by Zigrang-Sylvester correlation [10].The steady-state momentum balance of the closed loop can be expressed as The expressions for all the pressure drops in HEM model and drift flux model (except acceleration pressure drop in drift flux model) and the expression for average core void fraction in the HEM model are taken from van Bragt and van der Hagen [7].The expression for acceleration pressure drop in the drift flux model is taken from Kok and van der Hagen [4]: The above set of equations is converted into explicit boiling boundary equation and momentum equations ( 2) and ( 9).These two equations are solved using the FSOLVE function in MATLAB.

Formulation for stability
The governing equations for stability analysis are derived from one-dimensional unsteady-state equations of mass, momentum, and energy balances (1).These equations are derived using HEM model as discussed in [7,11] .The core is divided into single-phase and two-phase regions as explained in Section 3.1.Integrating the energy equation from core inlet to boiling boundary (Z bb ), the state equation for the boiling boundary is By integrating the continuity equation, we get the state equation of core average void dynamics given by The relation between core inlet mass flux density and the exit mass flux density as given below is obtained by integrating the energy equation in the two-phase region: Assuming linear variation of void fraction along the core, the relation between the core average void fraction and the core exit quality is given by Riser is an adiabatic section, divided into N R nodes, of equal length.Integrating the differential continuity equation along the riser gives the riser void dynamics equation as Integrating the energy equation along each node in riser gives the relation between the mass flux densities of adjacent nodes: The dynamic pressure is continuous along the closed natural circulation loop and can be expressed as Each term in (17) is the summation of inertial, gravitational, frictional, acceleration, and minor losses terms.The detailed expressions of these terms are given in [7] The total dynamic pressure drop (summation of inertial, gravitational, acceleration, friction) can be evaluated by integrating the differential momentum equation along each component of the boiling system considered.The friction factor "f " in core and riser is evaluated using the Darcy-Weisbach friction factor formula.It is assumed that the feed water flow rate is equal to the steam flow rate to make the drum dynamics steady.The core top is considered as a single node, and the riser is divided into two nodes of equal length.The resulting mathematical model is a set of ordinary differential equations with state variables as the G ci , z bb , a c , α r1 , α r2 .The system is linearized at the fixed point (obtained from steady-state analysis), and the Jacobian matrix is evaluated numerically.The nature of the eigenvalues determines the stability of the system.If all the eigenvalues lie in the left half of the complex plane then the system is stable.

Effects of operating and geometric parameters
Numerical simulations of NCBWR were performed by the lumped parameter model as well as the RELAP5 model, and the results were compared.The effects of various operating parameters like pressure, power, feedwater inlet enthalpy (FIE), and geometric parameters like riser length, core length, hydraulic diameter of core, etc were studied.While performing the parametric studies, certain parameters were varied, and other parameters were fixed.The input data are listed in Tables 1 and 2. The input parameters are chosen based on the information available in the literature [6] and suitable assumptions.The trends obtained are discussed in the following paragraphs.
At constant pressure and FIE, as power increases, the mass flow rate decreases (see Figure 3), except for the low power range, in which the mass flow rate increases with power.As power increases, the steam quality should also increase.Due to the presence of more vapor, two effects will occur.Firstly, with an increase in the steam quality, the density difference between the riser section and the downcomer section will increase.This, in turn, will tend to increase the mass flow rate.Secondly, the velocity of the coolant will increase as the steam quality increases.This will increase the frictional resistance and hence tends to reduce the mass flow rate.In most of the power range considered in the present study, the frictional resistance seems to be more dominant, so the mass flow rate decreases as shown in Figure 3.
At constant power and FIE, with an increase in the riser length, the mass flow rate increases (see Figure 4).As the riser length increases, two effects will occur.Due to longer riser length, the frictional resistance to the flow increases that tends to decrease the mass flow rate in the system.However, there is a gain in the buoyancy force due to the increased elevation.This will increase the flow rate in the system.The latter effect seems to be dominant.Hence, there is a notable rise in the mass flow rate of the system as the length of the riser increases.As the pressure and total power are kept constant, the mass flow rate decreases with an increase in core length (see Figure 5).As the core length increases, the mass flow rate may increase due to increased elevation The mass flow rate may decrease due to an increase in frictional resistance because of longer core.The flow area of the core is very small compared to the riser.Because of this, the frictional resistance is dominant; hence the mass flow rate decreases with an increase in core length.As the hydraulic diameter of the core increases, the frictional resistance should decrease due to an increase in the flow area, decrease in the aspect ratio, and decrease in the friction coefficient (due to the decrease in relative roughness).Therefore, the mass flow rate increases with hydraulic diameter of the core (see Figure 6).From Figures 3-6, it is observed that the values of mass flow rate predicted by the drift flux model are closer to those obtained with RELAP compared to the homogeneous equilibrium model.Therefore, drift flux model is the more accurate of the two lumped parameter models.For further analysis, lumped parameter model using the drift flux formulation is used.This is denoted as "LPM" in this paper.
As power and pressure are kept constant and FIE increases, the core inlet enthalpy (CIE) should also increase, as shown in Figure 7. Due to this, more sensible heat is available to the coolant, so the vapor generation rate increases.Consequently, the steam quality in the core and the riser should also increase and hence, the density of the mixture in the two-phase region should decrease.This, in turn, increases the frictional resistance in the two-phase region.The reduction in density will also increase the driving pressure difference due to buoyancy effects, which may result in an increased mass flow rate.The frictional resistance seems to be dominant.Therefore, as FIE increases, the mass flow rate decreases (see Figure 8).As power and pressure are kept constant and FIE increases, the void fraction should also increase as shown in Figure 9.This is because of the decrease in mass flow rate as the FIE is increased (see Figure 8).This will increase the heat removed per unit mass of the coolant.Due to this, the steam  quality increases; hence the void fraction also increases.At the same FIE and power, the void fraction decreases with an increase in pressure (see Figure 9).

Variations in critical heat flux (CHF)
CHF occurs during the transition from nucleate boiling to film boiling.Typically, this transition is accompanied G. V. Durga Prasad et al.  by a large temperature jump on the surface at which boiling occurs.This drastic increase in temperature can cause burnout or failure of equipment.Therefore, it is very important to know about the CHF, and when it will occur.In the present study, CHF is studied at the exit of the core.The chances for the occurrence of CHF are more at the exit of the core where the steam quality is high (at constant linear power).The presence of more vapor will impede the liquid flow back to the heated surface, and the liquid cannot rewet the heated surface.In the present analysis, the PG-CHF correlation [12] in the basic form is used.This basic form uses the local equilibrium quality and the local heat flux.In numerical simulations, the PG-CHF correlation is used as an alternative for the "CHF table look-up" method.This correlation uses the direct substitution method, which makes use of the available thermal-hydraulic and geometry information, and predicts the critical heat flux ratio (CHFR) at each point along the channel based on input heat flux at that point.
If the mass flux is high, the chances for the coolant to wet the heater wall are more.So, the transition from nucleate boiling to film boiling will occur at a relatively higher heat flux.If the steam quality is high, the vapor will surround the heater wall.So, the chances for coolant to wet the heater wall will be less.Therefore, the transition from nucleate boiling to film boiling will occur at a relatively lower heat flux.
In this analysis, the axial power profile is taken as uniform along the length of the core.From Figure 10, at fixed inlet conditions, with an increase in the power, CHFR (which is equal to the ratio of the critical heat flux to local heat flux) decreases.This is because, as power increases, the core exit steam quality increases (see Figure 10(a)), and the mass flux decreases (see Figure 10(b)), hence the transition from nucleate boiling to film boiling occurs quickly.So, the CHFR decreases with an increase in the power (see Figure 10(c)).As shown in Figure 11, with an increase in the core inlet enthalpy (CIE), the critical heat flux value first increases and then decreases.This is because, as CIE increases, the mass flux first increases and then decreases, and steam quality increases monotonically (see Figures 11(a) and 11(b)).At lower core inlet enthalpy, the mass flux increases, and steam quality also increases.But an increase in mass flux has the more influence on CHF than the increase in steam quality.So, the CHF value increases at lower core inlet enthalpy (see Figure 11(c)).At higher CIE, the CHF decreases with an increase in CIE.This is because, at higher CIE, mass flux decreases, and steam quality increases with an increase in CIE.Both effects contribute towards lower CHF as CIE increase.
From Figure 12(b), the CHF decreases with the core length.As the core length increases, the mass flow rate decreases (see Figure 5), and steam quality value increases (see Figure 12(a)).Because of both of these effects, the transition from nucleate boiling to film boiling becomes faster.So, CHF will occur at lower heat flux (see Figure 12(b)).With an increase in the hydraulic diameter of the active core, the CHF will increase as shown in Figure 12(b).As the hydraulic diameter increases, the mass flow rate also increases (see Figure 6).The quality decreases with an increase in the hydraulic diameter (see Figure 12(a)).Variation of both of these parameters will increase the CHF (see Figure 12(b)).From Figure 13, as the riser length increases, the CHF also increases.This is because, as the length of the riser increases, the mass flux increases (see Figure 4), and steam quality decreases (see Figure 13(a)).This will delay the transition from nucleate boiling to film boiling.So, the CHF will occur at relatively higher heat flux (see Figure 13(b)).

OPTIMIZATION OF PHT LOOP
Optimizing the material required to construct the PHT loop is very important from the economic point of view.In the present study, sequential quadratic programming (SQP) is adopted for the optimization procedure [13].The SQP algorithm has been one of the most successful general methods for solving nonlinear constrained optimization problems.Broadly defined, the SQP method is a procedure that generates iterates converging to a solution of the problem by solving a quadratic problem that is an approximation to the nonlinear problem.This type of method has been shown to be a very useful tool for solving nonlinear problems, especially where a significant degree of nonlinearity is present.

Formulation for optimization
The optimization procedure adopted here was to minimize the volume of the material required to construct the PHT loop of NCBWR.The design parameters for optimization technique are chosen as geometric parameters like length of the core, diameter of the coolant channel, length of the riser, diameter of the riser, and the cross sectional area of the downcomer.Appropriate upper and lower bounds were imposed on these design variables.The steady-state lumped parameter model ( 2)-( 9) was used for the calculations in the optimization procedure.
The heat flux to the coolant in the active core region should be less than the allowable heat flux (AHF) based on CHF considerations.If the applied heat flux crosses AHF, a stable vapor film may form on the surface of the core.If the heat flux to the core increases further, the fuel rod temperature will become very high.So, considering the critical heat flux constraint for optimization is important from safety point of view.
The system should be stable for the given operating conditions and at the minimum material condition.So, the dynamic stability of the system should be checked at the minimum material condition.The riser length is the most influencing parameter for natural circulation mass flow rate (see Figure 4).At the minimum material condition, the optimization code will adjust the length of the riser to the minimum possible value.At this situation, the mass flow rate of the system is very low.At this lower mass flow rate, the critical heat flux constraint will be satisfied, but the system may be unstable.To avoid the possibility of instabilities, a novel approach adopted here is to impose a minimum mass flux constraint in the optimization procedure.The mass flux in the active core region should be greater than or equal to a certain allowable mass flux (AMF) value.The value of AMF has to be chosen by trial and error, depending on the reactor power and operating conditions.If the resulting minimum material condition turns out to be unstable, AMF can be increased and another minimum material condition can be obtained.
Considering a typical steel material, the formula for the thickness of wall of the heater is taken as [14] T = (P + 17.26)D 1240 + 0.99 1000.
Here, pressure is in MPa, diameter is in mm, and thickness in meters.The volume of the material required for a pipe is Here, lengths, diameters, and thicknesses are in meters, volume is in m 3 .The optimization problem is where V c is the total volume of the core, V r is the total volume of the riser, V dc is the total downcomer volume.
Here, AHF is the allowable heat flux AHF = CHF/FOS, (FOS means the factor of safety), and X is a vector of design variables (geometric parameters).X min and X max are lower and upper bounds imposed on the design variables.AMF is the minimum allowable mass flux based on stability considerations.

Optimization results
Figure 14 shows the plot of minimum volume of the material required to fabricate the PHT loop versus power, at FOS of 1.4, and mass flux at 1000 kg/m 2 -s.It is seen that generally the volume of the material required increases with power.The variation of the material volume with the power is almost linear.The real part of the least stable eigenvalue of the Jacobian matrix has been indicated for minimum material condition at different powers in Figure 14.Since the least stable eigenmode lies in the left half plane (LHP), the system is stable for the given range of powers, at the minimum material condition.Yan et al. [15] have studied the effect of mass flow rate on the instabilities of natural circulation boiling water reactors.The mass flow rate in natural circulation depends on the geometric parameters and power.Therefore, the power and geometric parameters were taken as main governing factors for stability analysis.From Figure 15, it can be seen that, at the power of 220 MW, and for the minimum material condition, the system is stable at high mass fluxes, and unstable for the mass flux below 400 kg/m 2 -s.From this figure, it can be inferred that the system is unstable at low mass fluxes, if designed for the minimum material condition.Therefore, the use of the minimum mass flux constraint in the optimization procedure is justified.The system has largest stability margin for the mass flux of 800 kg/m 2 -s.
In Figure 16, the material required versus factor of safety has been plotted at power of 220 MW and mass flux of 1000 kg/m 2 -s.At lower values of the FOS, the volume of the material required is almost constant.At higher values of FOS, the volume of the material required is very sensitive to FOS.The volume of the material required increases drastically at higher FOS.The system is stable for the range of FOS considered.
The factor of safety should be chosen to provide sufficient safety margin in the presence of uncertainty in CHF prediction and other model uncertainties If FOS is chosen as 1.5, then the resulting design will be reasonably safe from the point of view of CHF as well as stability margin.However, it will be less economical compared to the design based on a lower FOS, for example, 1.4.
Table 3 shows parameters corresponding to 3 designs.Design 1 corresponds to FOS of 1.4 and mass flux of 1000 kg/m 2 -s.Design 2 corresponds to FOS of 1.4 and mass flux of 800 kg/m 2 -s.Design 3 corresponds to FOS of 1.5 and mass flux of 1000 kg/m 2 -s.It is seen that design 2 and design 3 are better than design 1, based on stability margin criterion (see Table 3).Out of designs 2 and 3, from the point of view  of economy, design 2 is better, whilst from the point of view of safety and stability margin, design 3 is better.Thus, the optimization problem requires higher level decisions based on Pareto-optimality [16].

CONCLUSIONS
Two lumped parameter steady-state models of NCBWR have been developed based on integral balance equations of mass, momentum, and energy, using the drift flux model and the HEM model for two-phase flow.Numerical simulations have been carried out with these models as well as RELAP5/MOD3.4code.The parametric trends obtained from the lumped parameter models are found to be similar to those calculated with RELAP5.Of the two lumped models, the one based on the drift flux model predicts values closer to those computed with RELAP5.The variation of CHF with different geometric and operating conditions is investigated.The following parametric trends are observed.
(i) With an increase in power, the steam quality increases, the mass flux decreases, and CHFR decreases.
(ii) With an increase in CIE, the steam quality increases, the mass flux first increases and then decreases, and CHF also first increases and then decreases.
(iii) With an increase in the core length, the steam quality increases and the mass flow rate decreases, hence CHF decreases.
(iv) With an increase in the hydraulic diameter of the core, the mass flow rate increases and the steam quality decreases, hence CHF increases.
(v) With an increase in the riser length, the mass flow rate increases and the steam quality decreases, hence CHF increases.
The optimum design of the PHT loop is obtained using the SQP algorithm based on safety and stability constraints.

Figure 3 :
Figure 3: Variation of mass flow rate with power.

Figure 4 :Figure 5 :
Figure 4: Variation of mass flow rate with riser length.

Figure 6 :
Figure 6: Variation of mass flow rate with hydraulic diameter of the core.

Figure 8 :
Figure 8: Variation of mass flow rate with FIE.

Figure 9 :
Figure 9: Variation of void fraction with FIE.

Figure 10 :
Figure 10: The variations of (a) steam quality, (b) mass flux, and (c) CHFR with power.

Figure 11 :
Figure 11: The variations of (a) steam quality, (b) mass flux, and (c) CHF with CIE.

Figure 12 :
Figure 12: The variations of (a) steam quality and (b) CHF with core length.

Figure 13 :
Figure 13: The variations of (a) steam quality and (b) CHF with core length.

Figure 14 :
Figure 14: Variation of minimum volume of the material required with power.

Figure 15 :
Figure 15: Variation of minimum volume of the material required with mass flux.

6 Figure 16 :
Figure 16: Variation of minimum volume of the material required with FOS.

Table 3 :
Design parameters at minimum material condition.