Development of an improved non-equilibrium multi-region model for pressurized water reactor pressurizer

Annals of Nuclear Energy 126 (2019) 133–141

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Development of an improved non-equilibrium multi-region model for pressurized water reactor pressurizer Xianping Zhong a,⇑, Xiaolong Zhang a, Jiyang Yu a, Muhammad Saeed a, Yi Li b, Zhihui Chen b, Bin Tang b, Yan Sun b a b

Department of Engineering Physics, Tsinghua University, Beijing 100084, China Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu 610213, China

a r t i c l e

i n f o

Article history: Received 1 August 2018 Received in revised form 10 September 2018 Accepted 4 November 2018

Keywords: Non-equilibrium multi-region model Pressurizer Gauss-Jordan elimination with partial pivoting Mesh-independent solution

a b s t r a c t This paper presents an improved non-equilibrium multi-region model for accurate prediction of pressure in the pressurizer of pressurized water reactor under transient conditions. The pressurizer model consists of three layers: a liquid layer, a saturated layer and a vapor layer. Each layer is further meshed into several control volumes (regions). The mathematical model derived from the mass and energy conservations describes most of the important thermal-hydraulics phenomena that occur in the pressurizer: bulk evaporation (bubble rise), bulk condensation (condensate fall) and spray condensation. The surge flow caused by the temperature change of the primary coolant is treated as the external boundary condition. To obtain the numerical solutions, the control equations of this model are first transformed into a set of linear equations, which is subsequently solved by the Gauss-Jordan elimination with partial pivoting. The mesh-independent solutions of the pressure of this model are in good agreement with the Shippingport pressurizer loss-of-load experiments. Moreover, the transient flow field of the pressurizer is analyzed. Unlike most of the previous models that require changing the empirically-determined surge flow distribution coefficient(s) to match the experiment data, this model does not need the knowledge of the surge flow distribution. And unlike the previous non-equilibrium multi-region model, this model can obtain the mesh-independent solutions. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction The pressurizer is one of the most essential components in the Pressurized Water Reactor (PWR), that maintains the pressure of primary coolant system in steady-state operation, limits the pressure change within permissible value during operational transient and prevents the primary coolant system from overpressure in accidental scenarios. Under normal conditions, the pressurizer consists of a lower liquid lump and an upper vapor lump. The lower liquid lump is usually connected to the hot leg of the primary circuit through a surge line and the upper vapor lump to the cold leg of the primary circuit through a spray line. Temperature changes in the primary coolant can result in water surges in or out of the lower liquid lump. During an insurge, spray may be activated to condense steam to provide room for entering water and the electric heater be activated to maintain the water in saturated. During an outsurge, the electric heater may be activated to prevent pressure from decreasing owing to the departure of the ⇑ Corresponding author. E-mail address: [email protected] (X. Zhong). https://doi.org/10.1016/j.anucene.2018.11.010 0306-4549/Ó 2018 Elsevier Ltd. All rights reserved.

water and the sudden expansion of the steam. With the drastic pressure changes in the pressurizer, the liquid and vapor lumps may turn into two phase state. As a result, the bulk evaporation and bulk condensation may occur in the pressurizer, which in turn affects the pressure of the pressurizer. Since the pressure of the pressurizer is essential for the reactor performance, it is required to predict the pressure response of the pressurizer with great accuracy during transients. The pressurizer dynamic analysis models, depending on whether the separated liquid and vapor phases in the pressurizer are assumed to be at the saturation state, are classified into two categories: equilibrium models and non-equilibrium models (Baggoura and Martin, 1983). The equilibrium models assume that the liquid and vapor phases in the pressurizer are saturated at all times. The governing equations of these models are quite simple and the calculation speed is fast. Some pressurizer programs based on these models have been applied to the simulation of the overall nuclear plant system (Christy and Galan, 1973; Thakkar, 1975). The equilibrium models are an approximation to the actual working condition. When the pressure changes, the liquid and vapor phases of these

134

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

Nomenclature A B C c1 c2 cCR cp db g h
h hsu hfg hsp Ja LLa m n P QH QTC T t VBR VCF VL

cross section area of the pressurizer [m2] sparse coefficient matrix of Eq. (28) constant term vector of Eq. (28) constant of Eq. (14) constant of Eq. (14) constant of Eq. (16) specific heat at constant pressure [J/(kg*K)] bubbles rising diameter [m] gravity acceleration [m/s2] specific enthalpy [J/kg] flow specific enthalpy [J/kg] specific enthalpy of insurge flow [J/kg] latent heat of vaporization of water [J/kg] inlet specific enthalpy of the spray [J/kg] Jakob number Laplace length [m] number of the control volumes in liquid layer number of the control volumes in vapor layer pressurizer pressure [Pa] heat absorption rate from the heater [W] Fourier thermal conductivity absorption rate from adjacent control volume(s) [W] temperature [K] time [s] bubble rise velocity [m/s] condensate fall velocity [m/s] pressurizer liquid volume [m3]

models will deviate from their original state and gradually reach a new equilibrium state through the heat and mass exchange between each other phase. The predictions by such models are satisfying as long as the transient processes are slow. However, these equilibrium models do not meet the requirement for accurate predictions under rapid transients, such as the small break loss of coolant accident and the load rejection accident of steam turbine. The non-equilibrium models consider the vapor and liquid phases by applying conservation laws (mass, energy and perhaps momentum) on each control volume of the pressurizer. The vapor phase can be in superheated or two-phase saturated state, and the liquid phase in subcooled or two-phase saturated state. The mass and energy exchange between two phases can be achieved through condensate fall and bubble rise. According to the number of the control volumes employed for the pressurizer, the nonequilibrium models are further classified as two-region models, three-region models and multi-region models. In other words, a region corresponds to a control volume. The non-equilibrium two-region models employ only two control volumes for each phase (Abdallah et al., 1982; Baron, 1973; Kim et al., 2006; Nahavandi and Makkenchery, 1970; Redfield et al., 1968). Most of the above models assume that the insurge flow mixes thoroughly with the liquid volume (see Fig. 1A). The nonequilibrium three-region models are considered as an improvement of most of the non-equilibrium two-region models, which employ two liquid volumes for the liquid phase because of the stratification caused by the insurge subcooled coolant. The non-equilibrium three-region model developed by Baggoura and Martin (1983) employs an upper water volume and a lower water volume empirically. They assume that the surge flow mass rate is dispensed into these two water volumes simulataneously in the ratio expressed by an empirical distriubtion factor b (see Fig. 1B). The nonequilibrium three-region model developed by Min Baek et al. (1986) employs a upper saturated water volume and a lower surge

VT W WBC WBC sum WBE WBEsum WSC WSC sum W rsv W sp W su X

a

b

ch cp q l r t

v fg x

pressurizer total volume [m3] mass flow rate [kg/s] bulk condensation mass rate [kg/s] total bulk condensation mass rate from the vapor layer [kg/s] bulk evaporation mass rate [kg/s] total bulk evaporation mass rate from the liquid layer [kg/s] vapor condensation mass rate by the spray droplets [kg/s] total vapor condensation mass rate by the spray droplets [kg/s] mass flow rate of relief and safety valve flow [kg/s] mass rate of spray water [kg/s] mass flow rate of surge flow [kg/s] solution vector of Eq. (28) void fraction empirical distribution factor for surge flow partial derivative of Eq. (24) [m3/J] partial derivative of Eq. (24) [m3/(kg*Pa)] density [kg/m3] dynamic viscosity [Pa*s] surface tension [N/m] specific volume [m3/kg] difference between the saturated steam specific volume and saturated water specific volume [m3/kg] thermodynamic quality of the control volume

water volume empirically. The lower water volume absorbes the whole insurge water or ejectes the whole outsurge water (see Fig. 1C). The non-equilibrium multi-region model is originally developed by You et al. (2001). The model divides the pressurizer into a liquid layer, a saturated layer and a vapor layer vertically. The liquid layers and the vapor layers consist of several control volumes. The saturated layer is formed of a saturated water volume and a saturated steam volume near the vapor-liquid interface. By doing this, the model can be used to analyze the flow field of the pressurizer. However, the model also presumes that the surge flow mass rate is distributed into each liquid volume simulataneously and the distribution factors bi are determined empirically (see Fig. 1D). In addition, the simulation results are mesh-dependent (You et al., 2001). Although the non-equilibrium three-region and multi-region models mentioned above can predict the pressurizer performance well, most of them require changing the empirically-determined surge flow distribution coefficient/coefficients to match the experiment data, which brings/bring great uncertainty to the models. This study assesses the characteristics of the aforementioned models, develops an improved non-equilibrium multi-region model to predict the pressure response of the pressurizer during transients and verifies this model by comparing the Shippingport pressurizer loss-of-load experiments. There are two novelties of this model. Firstly, this model can simulate the propagation of the surge flow from the surge line inlet to the upper water layer. Therefore, there is no need of the knowledge of the surge flow distribution. Secondly, this model can obtain the mesh-independent solutions, which is an improvement of the previous nonequilibrium multi-region model. 2. Model development The pressurizer consists of three layers: a liquid layer, a saturated layer and a vapor layer. It is divided into m liquid control

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

135

Fig. 1. Schematics of non-equilibrium two-region model (A), non-equilibrium three-region model developed by Baggoura and Martin (B) and by Min Baek et al. (C), nonequilibrium multi-region model developed by You et al. (D), and the present model (E).

volumes and n vapor control volumes (see Fig. 1E). The liquid control volumes 1 to m-1 belong to the liquid layer. The vapor control volumes 1 to n1 belong to the vapor layer. The m’th liquid control volume and the n’th vapor volume near the liquid-vapor interface are considered to be in saturated state. These two volumes form the saturated layer. Assumptions that are made in the model include:

7. Delay times of bubbles rise and condensate fall are neglected. 8. Wall condensation and boiling are neglected. 9. Work and gravitational potential energy of each control volume are neglected. 10. Positive direction of the liquid control volume is upward and positive direction of the vapor control volume is downward.

1. Each control volume has the same averaged uniform cross section area. 2. Each control volume has the same pressure. 3. Each control volume has the averaged uniform specific enthalpy at each time step. 4. Both the liquid spray and the condensation on the liquid spray become saturated before reaching the liquid-vapor interface. 5. Transient time for liquid spray droplets reaching the liquidvapor interface is neglected. 6. Mass exchange at the liquid-vapor interface is a result of bulk evaporation (bubble rise) and bulk condensation (condensate fall).

In this work, physical variables of the control volumes consist of state variables and flow variables. The state variables are variables that represent the properties of the control volume itself, such as density, specific enthalpy, bulk evaporation mass rate and bulk condensation mass rate. The flow variables are variables that represent the transport values of the two adjacent control volumes. The mass flow rate (denoted by W) and the flow specific enthalpy


(denoted by h) are the only two flow variables in this model. By X l;i we denote the state variable X in the i’th liquid control volume or the flow variable X between the i’th liquid control volume and the (i + 1)’th liquid control volume. By X v;i we denote

136

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

the state variable X in the i’th vapor control volume or the flow variable X between the i’th vapor control volume and the (i + 1)’th vapor control volume. By X sf and X sg we denote the physical variable X at the saturated water state and the saturated steam state, respectively.



 L d ðqsf hsf  PÞ VmL þ qsg hsg  P V T V n dt





¼ WBEsum hsg þ WBC sum hsf þ W l;m1 h l;m1 þ W v;n1 h   þ W sp þ WSC sum hsf þ QTC l;m þ QTC v;n :

v;n1

ð12Þ

2.1. Governing equations 2.2. Physical models In the liquid layer, the time-dependent mass and energy conservation equations for volume i (i = 1,2,. . .,m-1) have the following forms:


 d ql;i VmL

¼ W l;i1  W l;i  WBEl;i ; ð1Þ dt  

d ql;i hl;i  P VmL

¼ W l;i1 hl;i1  W l;i hl;i  WBEl;i hsg þ QHl;i þ QTC l;i : dt ð2Þ The surge line is considered as a virtual liquid control volume denoted by ‘‘0”, thus the mass flow rate between liquid control volume 0 and liquid control volume 1 is given by:

W l;0 ¼ W su:

ð3Þ

By doing this, the surge flow (as a forcing function of the governing equations) can be spread to the upper liquid control volumes without introducing the surge flow distribution factor(s). The flow specific enthalpy between the liquid control volumes is determined by the donor control volume and thus it is given by:






hl;0 ¼ 




hl;i ¼

hsu

for W l;0 > 0

hl;1

for W l;0 < 0

;

hl;i

for W l;i > 0

hl;iþ1

for W l;i < 0

ð4Þ

WBEl;i ¼ qsg al;i A  VBRl;i ; ði ¼ 1; 2; . . . ; m  1Þ :

ð5Þ

In the vapor layer, the mass and energy conservation equations for volume i (i = 1,2,. . .,n1) are given as:


 L d qv;i V T V n dt

¼ W v;i1  W v;i  WBC v;i  WSC v;i ;



 L d qv;i hv;i  P V T V n dt




ð6Þ




¼ W v;i1 hv;i1  W v;i hv;i  WBC v;i hsf  WSC v;i hsg þ QTC v;i :

ð7Þ

Similarly, the relief and safety valve is considered as a virtual vapor control volume denoted by ‘‘0”, thus the mass flow rate between vapor control volume 0 and vapor control volume 1 is given by:

W v;0 ¼ W rsv :

ð8Þ

The flow specific enthalpy between the vapor control volumes is determined by the donor control volume and thus it is given by:


hv;0 ¼ hv;1


hv;i ¼



ð*W v;0 ¼ W rsv

hv;i

for W v;i > 0

hv;iþ1

for W v;i < 0

 0Þ ; ði ¼ 1; 2; . . . ; n  1Þ :

ð10Þ

dt

þ W sp þ WSC sum ;

8 gðqsf qsg Þd2b > > < c1 lsf  0:25 VBR ¼ gð q > sf qsg Þrsf > : c2 q2 sf

ð11Þ

for LLa < db for LLa  db

;

ð14Þ

where c1 is related to the physical properties of the liquid phase and between 2/9 and 3/9. c2 is experimentally determined to be 1.18. The Laplace length LLa is determined by:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

LLa ¼

rsf

gðqsf  qsg Þ

:

ð15Þ

The bubbles rising diameter db is deemed to be equal to the bubbles departure diameter in boiling saturated liquids, which is given by Cole and Rohsenow correlation (Cole and Rohsenow, 1969):

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

db ¼ cCR Ja

rsf

1:25

gðqsf  qsg Þ

;

ð16Þ

where cCR is depended on the liquid type and is equal to 1.5E4 for water. The Jakob number Ja is expressed as:

Ja ¼

¼ WBEsum þ WBC sum þ W l;m1 þ W v;n1

ð13Þ

where the bubble rise velocity VBR is determined from Gunther and Kreith correlation (Gunther and Kreith, 1950):

ð9Þ

In the saturated layer, the mass and energy conservation equations for volume are given as:


 L d qsf VmL þ qsg V T V n

2.2.1. Bulk evaporation and condensation During the pressure change in the pressurizer, a part of the water in the liquid layer becomes superheated and flashes into the saturated layer and a part of the steam in the vapor layer becomes subcooled and condenses into the saturated layer. One method that describes these phenomena is applying the kinetic theory of gases and using experimentally determined evaporation and condensation coefficients (Baggoura and Martin, 1983). Another method to model the bulk evaporation and condensation is assuming that these phenomena come from bubble rise and condensate fall. Abdallah et al. (1982) suggest that the mass rate of the bubble rise and condensate fall is equal to the time rate of change of the bubble mass in the liquid layer and the liquid mass in the vapor layer, respectively. However, most of the pressurizer transient analysis models (Baron, 1973; Min Baek et al., 1986; Nahavandi and Makkenchery, 1970) assume that the bulk evaporation and condensation is dependent on the velocity of the bubble rise and condensate fall, respectively, which is also adopted in this work. The mass rate of bubbles escaping from the liquid control volume i (i = 1,2,. . .,m1) is determined by:

qsf cpsf T sf : qsg hfg

ð17Þ

As the delay times of bubbles rise are neglected, the mass rate of bubbles entering into the saturated layer is determined by:

WBEsum ¼

m1 X i¼1

WBEl;i :

ð18Þ

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

The mass rate of condensate escaping from the vapor control volume i (i = 1,2,. . .,n1) is determined by:

WBC v;i ¼ qsf ð1  av;i ÞA  VCF v;i ;

ð19Þ

where the condensate fall velocity VCF is given as a constant. Similarly, as the delay times of condensate fall are neglected, the mass rate of condensate droplets entering into the saturated layer is determined by:

WBC sum ¼

n1 X

WBC v;i :

ð20Þ

where

v fg ¼ v sg  v sf :

BX ¼ C;

ð28Þ

where the solution vector



The homogeneous-fluid model is applied in both the liquid and vapor control volume. Therefore, the void fraction of these control volumes is given by:

8 1 for 0 < x < 1 > q > < 1þqsgsf ð1x1Þ for x  1

1 > > : 0

:

ð21Þ

for x  0

2.2.2. Condensation on spray The spray-condensate mixture is considered as an addition to the saturated layer, since the assumption is made that the mixture becomes saturated before reaching the liquid-vapor interface without a time delay. The total mass rate of vapor condensation on the spray droplets is given by:

WSC sum ¼ W sp

ð27Þ

After the expansion and rearrangement of Eqs. (1), (2), (6), (7), (11), (12), one can derive a linear equation set of 2(m + n  1) equations:

i¼1

a¼

hsf  hsp : hfg

ð22Þ

dp dV L dhl;1 dhl;2 dhl;m1 dhv;1 dhv;2 dhv;n1 ; ; ; ;...; ; ; ;...; ; dt dt dt dt dt dt dt dt T W l;1 ; W l;2 ; . . . ; W l;m1 ; W v;1 ; W v;2 ; . . . ; W v;n1 ð29Þ

X¼

2 R2ðmþn1Þ1 , the coefficient matrix B 2 R2ðmþn1Þ2ðmþn1Þ , and the constant term vector C 2 R2ðmþn1Þ1 . Eq. (28) is solved by the Gauss-Jordan elimination with partial pivoting (Cheney and Kincaid, 2012). 2.4. Code flow chart A code named RPDCODE is developed based on the governing equations, physical models and numerical method discussed above. The water/steam properties calculated in RPDCODE are extracted from the IAPWS-IF97 (Wagner and Kretzschmar, 2008). The calculation steps of RPDCODE are outlined below and illustrated schematically in Fig. 2.

It is hard to predict the mass of condensate of every vapor control volumes because of the complicated energy and mass transfer of the falling spray droplets. This work assumes that each vapor control volume of the vapor layer contributes the same to the total mass of the vapor condensation, which is expressed as: sum WSC v;i ¼ WSC n1

ði ¼ 1; 2; . . . ; n  1Þ:

ð23Þ

2.2.3. Pressurizer components The pressurizer components include the electric heater, the spray and the safety and relief valve. These components in the Shippingport pressurizer of the following validation are simply on-off type (see Table 1). More complex control logics and performance characteristics of these components have been widely discussed (Abdallah et al., 1982; Baggoura and Martin, 1983; Baron, 1973; Min Baek et al., 1986; Nahavandi and Makkenchery, 1970), which are essential to the model of the pressurizer but not repeated in this work. 2.3. Numerical method The time derivative terms in Eqs. (1), (2), (6), (7), (11), (12) are expanded into total derivative with respect to the time term t. In particular, the time derivative of density is expanded into:



  dq @ q

dh @ q

dp 1 dh dp ; þ

¼  2 ch þ cp ¼

dt dt dt @h p dt @p h dt v

ch ¼

( @v

@h

v fg hfg

137

for single phase for two phase

;

8 for single phase < @@pv ; cp ¼
dv fg dv sf  v fg
dhfg dhsf  : x þ dp  h x dp þ dp for two phase dp

ð24Þ

ð25Þ

ð26Þ

fg

Fig. 2. Code flow chart of RPDCODE.

138

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

1. Input the pressurizer geometry parameters (e.g. the pressurizer height and cross section area), model empirical coefficients (e.g. the condensate fall velocity), boundary conditions (e.g. control logics of the surge line, the heater, the spray and the relief and safety valve) and initial conditions (e.g. the pressurizer pressure and water level). 2. Update the state of the pressurizer control volumes based on the input data at the first simulation time and based on the solution X of Eq. (28) in the remaining simulation time. 3. Update the boundary conditions, namely the responses of the surge line, the heater, the spray and the relief and safety valve. The boundary conditions are updated by the external and internal factors. The surge flow rate is the external factor treated as an external forcing function while the responses of the heater, the spray and the valve are the internal factors controlled by the pressurizer state (e.g. the pressurizer height and water level). 4. Update the coefficient matrix B and constant term vector C of Eq. (28) based on the state of the control volumes and the boundary conditions. 5. Solve Eq. (28) by the Gauss-Jordan elimination with partial pivoting to obtain the solution vector X. 6. Make a judgment whether the simulation time is over. If yes, finish the code. If no, update the simulation time and go back to step 2.

Fig. 3. Calculation and experiment of pressure of the pressurizer due to loss-of-load from 74 MW(e).

3. Validation 3.1. Experiment introduction RPDCODE is verified with two Shippingport pressurizer loss-ofload experiments. The pressurizer is cylindrical and the pressurizer is maintained by three bands of electrical heaters, a spray and a steam relief valve. There are two pressure detectors used in the experiments. One is a narrow-range detector that reads from 12.06 to 15.51 MPaG and has an absolute pressure uncertainty of 3.45E2 MPa. The other is a wide-range detector reads from 0 to 20.68 MPaG and has an absolute pressure uncertainty of 0.21 MPa. One test is initiated from a reactor power level of 74 MW(e), reduced to 10 MW(e) and then to zero. The other is initiated from a reactor power level of 105 MW(e), reduced to 35 MW (e), 10 MW(e) and then to zero. The power imbalances of these two tests resulted in the temperature changes of the primary coolant, which further induce the combination of the insurge, outsurge, spray actuation and heater cycling in the pressurizer (Redfield et al., 1968). The input data for the code simulations are listed in Table 1 (Nahavandi and Makkenchery, 1970; Redfield et al., 1968; Sun,

Fig. 4. Calculation and experiment of pressure of the pressurizer due to loss-of-load from 105 MW(e).

2004). Two of the three banks of the heaters are not triggered in the experiments and are not triggered in the following simulations. The steam relief valve is even not applied in the experiments.

Table 1 Input data for the code simulations. Parameter type

Parameter

Geometry parameters Model empirical coefficients Boundary conditions

Heater

Spray

Surge line Initial conditions

Value/source Inner diameter (m) Height (m) c1 of Eq. (14) (–) VCF of Eq. (19) (m/s) Thermal power (kW) Turn-on setpoint (MPa) Turn-off setpoint (MPa) Volume flow rate (m3/s) Inlet temperature (K) Turn-on setpoint (MPa) Turn-off setpoint (MPa) Mass flow rate (kg/s) Insurge flow temperature (K) Saturated pressure (MPa) Water level (m)

1.37 5.61 2/9 5.0 80.0 (Distributed equally to each control volume in the liquid layer) 13.655 13.793 2.397E-3 537.15 (Cold leg temperature) 14.276 13.966 Nahavandi and Makkenchery (1970) 555.15 (Hot leg temperature) Extracted from test curves 2.482

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

Therefore, the parameters of the above components are not shown in the Table 1. It should be noted that the time histories of the surge flow rate of the experiments, as the external factors of the boundary conditions, are not available to the authors. By analysis of the time histories of the water level and the spray flow rate, Nahavandi and Makkenchery (1970) have inferred the surge mass flow rate time histories of these two experiments for their calculation input, which are adopted but not repeated in this study. Similarly, since the authors are uncertain about the geometry and the location of the heaters, the thermal power of the heaters is assumed to be distributed equally to each of the control volumes in the liquid layer.

139

3.2. Results and discussion The comparisons between the simulations and the experiments of the pressure transients of the loss-of-load from 75 MW(e) and 105 MW(e) are shown in Figs. 3 and 4, respectively. Figs. 3 and 4 show that when the number of control volumes increases to 400 (m = 200, n = 200) or higher, the pressure behavior predicted by the code converges to a stable curve. The absolute errors of the pressure between the mesh-independent simulations and the experiments are less than 0.5 MPa. The transient flow fields of liquid and vapor layers are shown in Figs. 5–8. The number of control volumes is set to 500 (m = 250,

Fig. 5. Mass flow rate between liquid control volumes due to loss-of-load from 74 MW(e) (m = 250, n = 250).

Fig. 6. Mass flow rate between vapor control volumes due to loss-of-load from 74 MW(e) (m = 250, n = 250).

140

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

Fig. 7. Mass flow rate between liquid control volumes due to loss-of-load from 105 MW(e) (m = 250, n = 250).

Fig. 8. Mass flow rate between vapor control volumes due to loss-of-load from 105 MW(e) (m = 250, n = 250).

n = 250) to obtain the mesh-independent results. As can be seen from Figs. 5 and 7, the mass flow rate of the liquid layer decreases with the distance from the surge line inlet, which reflects the propagation of surge flow in the liquid layer without introducing the surge flow distribution factor(s). It is noteworthy that due to the bulk evaporation, the mass flow rate in the liquid layer sometimes fluctuates dramatically (see
80 s of Fig. 5 and
480 s of Fig. 7). The new model can yield mesh-independent results, which the preceding model created by You et al. (2001) fails to achieve. The model approximates the experimental data quite well. The distribution factors for surge flow that are indispensable for previously

mentioned models are not included. Since the factors are casedependent and cannot be derived from general rules, this gives the new model good versatility and expansibility. Two possible reasons are proposed to explain the discrepancy between the simulations and experimental data. The first is the dearth of firsthand input parameters for the code. e.g., the surge mass flow rate time histories of Shippingport experiments derived by Nahavandi and Makkenchery (1970) should be deemed approximate. The other cause of computational errors lies in the model of the bubble rise and condensate fall. The mechanism behind the phenomenon has not been thoroughly understood yet. Particularly, the

X. Zhong et al. / Annals of Nuclear Energy 126 (2019) 133–141

extreme conditions such as high pressure and rapid pressure changes, make it more cumbersome to predict the velocity of condensate fall precisely. As a result, the model introduces a uniform empirical-determined condensate fall velocity to solve the problem. Ideally, the most reliable and accurate prediction of the behavior of the bubble rise and condensate fall can be obtained by direct and numerical solution of the Navier-Stokes equations while it requires enormous computer calculation resources. The model in this work is an appropriate option for developing a fast response pressurizer program. 4. Conclusions This paper presents an improved non-equilibrium multi-region model. The mathematical model divides the pressurizer into a liquid layer, a saturated layer and a vapor layer, which are further meshed into several control volumes. The model derived from the mass and energy conservations describes most of the important thermal-hydraulics phenomena occurring in the pressurizer: bulk evaporation (bubble rise), bulk condensation (condensate fall) and spray condensation. The model introduces a virtual liquid control volume to simulate the surge line and treats the surge flow as the mass flow between two adjacent liquid control volumes. Therefore, the surge flow can be spread to upper liquid layer through the governing equation without needing the knowledge of the surge flow distribution. The model can obtain the mesh-independent solutions and the absolute errors of the pressure between the mesh-independent simulations and Shippingport pressurizer loss-of-load experiments are less than 0.5 MPa, which verifies the reliability of the model. Acknowledgement This project is supported by the Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China.

141

References Abdallah, A.M., Mariy, A.H., Rabie, M.A., Nagy, M.E.S., 1982. Pressurizer transients dynamic model. Nucl. Eng. Des. 73, 447–453. Baggoura, B., Martin, W.R., 1983. Transient analysis of the three Mile island unit 2 pressurizer system. Nucl. Technol. 62, 159–171. Baron, R.C., 1973. Digital model simulation of a nuclear pressurizer. Nucl. Sci. Eng. 52, 283–291. Cheney, E., Kincaid, D., 2012. Numerical Mathematics and Computing. Nelson Education. Christy, P.T., Galan, V.J., 1973. POWER TRAIN: General Hybrid Simulation for Reactor Coolant and Secondary System Transient Response (No. BAW–10070). Babcock and Wilcox Co.. Cole, R., Rohsenow, W.M., 1969. Correlation of bubble departure diameter for boiling of saturated liquids. Chem. Eng. Prog. Symp. Ser. 65, 211–213. https:// doi.org/10.1016/j.ijheatmasstransfer.2011.04.007. Gunther, F.C., Kreith, F., 1950. Progress Report. Jet Propuls. Lab. Calif. Inst. Technol.. Kim, T.W., Kim, J.W., Park, G.C., 2006. Development of nonequilibrium pressurizer model with noncondensable gas. Nucl. Eng. Des. 236, 375–384. https://doi.org/ 10.1016/j.nucengdes.2005.09.003. Min Baek, S., Cheon No, H., Park, I.Y., 1986. A nonequilibrium three-region model for transient analysis of pressurized water reactor pressurizer. Nucl. Technol. 74, 260–266. Nahavandi, A.N., Makkenchery, S., 1970. An improved pressurizer model with bubble rise and condensate drop dynamics. Nucl. Eng. Des. 12, 135–147. https://doi.org/10.1016/0029-5493(70)90002-6. Redfield, J.A., Prescop, V., Margolis, S.G., 1968. Pressurizer performance during lossof-load tests at Shippingport: analysis and test. Nucl. Appl. 4, 173–181. Sun, Z., 2004. Nuclear Power Equipment. Harbin Engineering University Press (in Chinese). Thakkar, J.G., 1975. Correlation of theory and experiment for the dynamics of a pressurized water reactor. Wagner, W., Kretzschmar, H.-J., 2008. IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. Int. Steam Tables Prop. Water Steam Based Ind. Formul. IAPWS-IF97, 7–150. You, H., Cui, Z., Cheng, Y., 2001. A nonequilibrium multi-region model for nuclear power pressurizer. Nucl. Power Eng. 22, 133–137 (in Chinese).