Development of a multi-compartment containment code for advanced PWR plant

Nuclear Engineering and Design 334 (2018) 75–89

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Development of a multi-compartment containment code for advanced PWR plant

T

⁎

Yongzheng Chen, Y.W. Wu , M.J. Wang, Y.P. Zhang, B. Tan, D.L. Zhang, W.X. Tian, S.Z. Qiu, G.H. Su State Key Laboratory of Multiphase Flow in Power Engineering, Department of Nuclear Science and Technology, Xi’an Jiaotong University, Xi’an City 710049, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Containment thermal-hydraulics Design basis accident Severe accident Buoyant plume Hydrogen distribution Film tracking Counter-current flow

HPR1000 (Hua-long Pressurized Reactor) Nuclear Power Plant is a type of third generation advanced PWR developed in China. In this study, a home-made multi-compartment containment code, ATHROC (Analysis of Thermal Hydraulic Response Of Containment), was developed to analyze the thermal-hydraulics and hydrogen behaviors in HPR1000 containment. Two types of bulk fluids, the atmosphere bulk fluid and the pool bulk fluid, were used to model the transient behaviors of multiphase flow. The atmosphere bulk fluid consists of steam, noncondensables and homogeneously dispersed liquid droplets. The code contains comprehensive models, including flow model, heat and mass transfer model, engineering safety feature model, etc. A plume model was implemented to assess buoyancy driven plume movement inside a compartment. A counter-current flow model was developed to evaluate the buoyancy driven bidirectional exchange flow through the junctions. Also, a film tracking model was developed to simulate film flow on the contiguous walls. Special methods such as using variable time steps and non-uniform nodalization for heat sinks were applied to speed up the computation. Code assessments were carried out by simulating several separated effects tests (Phebus FPT0 test, JAERI tests) and integral effects tests (CVTR tests, NUPEC M-7-1 test). Comparing results between code simulations and experimental data showed that the code was capable of providing reasonable predictions.

1. Introduction In the nuclear power plants, the containment system is the last barrier preventing fission products from being released to the environment. Much research has been done on the containment thermal hydraulics experiments and analytical assessments. In the early years, research mainly focused on the dynamic response of containments following a design basis LOCA (Loss Of Coolant Accident) or MSLB (Main Steam Line Break), which was important for the proper design of associated containment structures. Recently, hydrogen behavior during a core degrade accident has received considerable attention because it poses more serious threat to the containment integrity when no mitigations are available (Royl et al., 2000). International consensus is that a detailed knowledge of containment thermal hydraulics is necessary to analyze the effectiveness of hydrogen mitigation methods (Karwat et al., 1999). Though a large number of experiments concerning containment thermal hydraulics have been done, they are usually performed at only a reduced scale and under limited conditions. Conducting prototypical experiments in a full scale containment is not economically feasible. Thus, computer codes which involve various

⁎

Corresponding author. E-mail address: [email protected] (Y.W. Wu).

https://doi.org/10.1016/j.nucengdes.2018.05.001 Received 4 December 2017; Received in revised form 27 April 2018; Accepted 2 May 2018 Available online 11 May 2018 0029-5493/ © 2018 Elsevier B.V. All rights reserved.

phenomenological models have become an important tool to perform safety analysis. Lumped-parameter codes are widely applied to integrated containment analysis such as probabilistic risk assessments (PRAs) and PSAs because they involve various coupling phenomenological models and are fast-running. During the past few decades, various containment codes with different degrees of sophistication and experimental validations have been developed throughout the world. CONTEMPT-LT/ 028 was developed to predict the long-term behavior of PWR containment systems subjected to postulated LOCA conditions (Lin, 1979). The containment conditions are obtained by solving the mass and energy balance equations. Different physical processes were introduced at each time step. However, only one to four compartments can be modeled in the code. CONTEMPT4/MOD6 is the latest version of the CONTEMPT series (Lin et al., 1986). It provides the capability of multicompartments simulation to model a prototypical containment more realistically. In addition, an implicit algorithm was employed to improve the numerical stability of flow calculations. However, the droplet field was not considered in the two CONTEMPT versions. In the CONTAIN 2.0 code, vapor, liquid pool, droplets as well as special phases like

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Nomenclature

total internal energy of the gas in cell i, [J] change of gas internal energy in cell i, [J] specific internal energy of steam, [J/kg] specific internal energy of gas in the cell i, [J/kg] change of specific internal energy of gas in the cell i, [J/ kg] us,sat , u w,sat saturated specific internal energy of steam and water at steam partial pressure, [J/kg] Vg gas volume of the node, [m3] specific volume of steam and gas, [m3/kg] vs , vg steam generation rate from the pool surface, [kg/s] Ws Wrain rainout rate from the steam, [kg/s] gas flow rate at junction ij and at junction k, [kg/s] Wij , Wk Zij,i , Zij,j central level of junction ij measured from the bottom of cell i and cell j, [m] bottom level of the two junction ends relative to the floor Zi , Zj of cell i and cell j, [m] Z wi , Z wj water level in the cell i and cell j, [m]

Ui δUi us ui δui

Aij,0, Aij,n water-covered junction area on the donor side of the flow at the beginning of iteration and at current iteration, [m2] CD discharge coefficient of the junction, [−] constant volume specific heat capacity of non-condensable c vi gas i, [J/(kg·K)] constant pressure specific heat capacity of gas, and gas in cpg , cpj the cell j, [J/(kg·K)] d diameter of the droplets, [m] Ddiff diffusivity of the steam, [−] DH hydraulic diameter, [m] Eg , Ew emissivity of gas and wall surface, [−] g gravitational acceleration, [m2/s] specific enthalpy of the water in the pool, gas in the cell i, h w , hj [J/kg] specific enthalpy of the steam and condensation film, [J/ hs , hcon kg] h w,sat , hs,sat saturated specific enthalpy of the water and steam at pool pressure, [J/kg] mass transfer coefficient, [kg/(s·m2·Pa)] hm hrad the radiation heat transfer coefficient, [W/(m2·K)] hnc natural convective heat transfer coefficient, [W/(m2·K)] hc convection transfer coefficient, [W/(m2·K)] hcd condensation heat transfer, [W/(m2·K)] gas flow resistance coefficient at junction ij, [−] Kij L characteristic length of a surface, [m] mass of steam [kg] ms mi mass of non-condensable gas i, [kg] mass of water in the pool, [kg] mw md mass of the droplets, [kg] molecular weight of water, [g/mol] Mw moles of non-condensable gas i, [mol] ni P total pressure, [Pa] Psat (T ) saturated steam pressure at temperature T , [Pa] partial pressure of steam and non-condensable gases, [Pa] Ps , Pnc cell pressures measured at the bottom of cell i and cell j, Pi , Pj [Pa] R the general gas constant, [−] temperature of gas, [K] Tg gas temperature in cell i and cell j, [K] Ti , Tj temperature at the interface between gas and film, [K] Tf temperature of the wall surface, [K] Tw Td temperature of the droplets, [K] Ug total internal energy of the gas, [J]

Greek letters

δmi Δ θ λl , λ g μl , μg ρi , ρj ρwi , ρwj ρg , ρf

ρd ρs σ

change of gas mass in the cell i, [kg] difference of values, [-] inclined angle of the heat sink, [。] thermal conductivity of water and gas, [W/(m·K)] viscosity coefficient of water and gas, [Pa·s] density of gas measured at the bottom of cell i and cell j, [kg/m3] density of water in the cell i and cell j, [kg/m3] density of gas in the cell and gas at the wall surface, [kg/ m3] density of the droplets, [kg/m3] density of the steam, [kg/m3] stefan constant, [−]

Dimensionless Numbers

Gr Nu Pr Re Ra Sh Sc

aerosols, fission products, and core debris are considered to predict containment conditions during severe accidents or design basis accidents (Murata et al., 1997). However, it is unable to predict the detailed gas-mixing behaviors due to the lack of a detailed entrainment model for jets, plumes, and shear layers and also the lack of a molecular/ turbulent diffusion model for gases mixing across gas-gas interfaces. Moreover, the momentum convection within a compartment was neglected. COCOSYS code (Allelein et al., 2008; Klein-Heβling et al., 2000) is being developed to simulate DBAs and severe accident in the containments of light water reactors. The code includes three main modules that treat thermal hydraulic, fission products and corium behavior, respectively. These modules communicate with each other via parallel virtual machine. Preliminary application of COCOSYS code in industry has been done in recently years. However, code improvements and validation process are still ongoing. CAP code written in C++ (Hong et al., 2015) is being developed in Korea aiming to analyze containment thermal hydraulic behaviors. Hydrogen distribution can be analyzed but hydrogen combustion is not considered at present. It

Grashof number, [−] Nusselt number, [−] Prandtl number, [−] Reynolds number, [−] Rayleigh number, [−] Sherwood number, [−] Schmidt number, [−]

can run in standalone mode or linked mode with RCS analysis code. Preliminary code assessments showed that the code could provide reliable predictions in DBAs. Unfortunately, film tracking model was not included in the code. Moreover, further assessments of the code on forced convection are necessary because of its great importance in the initial stage of transients. ASTEC (Van Dorsselaere et al., 2010; Van Dorsselaere et al., 2005), MAAP (Liang et al., 2014; Sehgal, 2011), MELCOR (Gauntt et al., 2005) are similar integral codes that can simulate most aspects of in-vessel and ex-vessel severe accident phenomena. To simulate long term transients in acceptable computation time, these codes employ simple physical models and numerical methods. Therefore, simulation results of these codes have large uncertainties. GOTHIC code seems to be the most advanced integral containment code until now (Hong et al., 2015). Three primary fields (steam/gas mixture, continuous liquid, liquid droplet) and two secondary fields (mist, ice) were considered (George et al., 2001). Besides using the typical lumped-parameter method, it provides the capability of detailed intensive analysis by subdividing some compartments. 76

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

fluid within each cell, was used to accelerate the calculation. Mass and energy conservation equations for each cell, momentum equations for each junction, were solved to obtain the rates of mass and energy variations. Fig. 2 shows the flow chart of the code computation.

To meet the industrial requirements of designing Chinese HPR1000 power plants (Xing et al., 2016), a home-made containment code named ATHROC (Analysis of Thermal Hydraulic Response Of Containment) was developed. This code development is in the scope of a long term research project named “Longteng 2020” that starts from 2014 and is still ongoing. The project aims to develop a series of homemade computer codes that could simulate the unique design systems of HPR1000. The ATHROC code is mainly applied to predict containment thermal-hydraulics and hydrogen behavior. Analytical models available in the code include models for intercell flow, heat and mass transfer processes (e.g., convection, thermal radiation, conduction, condensation) and engineering safety features (e.g., containment sprays, fan coolers, hydrogen igniter and recombiner). Compared with existing containment codes, ATHROC includes an advanced film tracking model for passive containment cooling system in AP1000 NPPs. A buoyant plume model was developed to consider the inhomogeneity inside a compartment, thus reducing the simulation distortion resulted from the assumption of homogeneity in most lumped parameter codes. A counter-current flow model was added to the general junction flow model to consider buoyancy-driven exchange flow. Moreover, special designed method is applied to speed up the computation. For containment code validation, many experiments are available as described by Chin et al. (2014). In this study, the code was assessed by separated effects tests and integral effects tests, such as Phebus FPT0 test (Clément et al., 2003; Gyenes and Ammirabile, 2011; Layly et al., 1996; Müller and Toth, 2007), JAERI tests (Motoki et al., 1983), CVTR tests (Schmitt et al., 1970; Tills et al., 2008) and NUPEC M-7-1 test (Stamps, 1995).

2.2. Methods to speed up the computation Several specially designed methods were included to speed up the calculation, such as the application of quasi-steady state model, tabularized results or correlations prior to detailed models, non-uniform nodalization for heat sinks, variable time step according to the intensity of transient process. The method of non-uniform nodalization for heat sinks and variable time step is briefly presented in this section. Total number of nodes for a heat sink is limited to 40. A different nodalization scheme is used for steel and concrete heat sinks. Steel is divided into uniform node sizes and concrete into non-uniform node sizes. For concrete heat sinks, ten fine boundary nodes are allocated to each face of a two-sided wall. Since the thermal resistance of concrete is very large, the thickness of each fine node is set as a small value of 0.0025 m. Accordingly, the thickness of coarse nodes in the central part of the concrete could be easily calculated. This nodalization method ensures that the steep temperature gradient near the wall surface could be accurately traced when the concrete wall is subjected to a high heat flux. In the input file, users can specify the maximum allowed fractional change of key state variables in one time-step. ATHROC first calculates the rates-of-change of the state variables based on preceding time step Δt0 and preceding values of the state variables. Then, for each specified variable Xi , the maximum allowed fractional change FXi and preliminary calculated rates-of change Xi̇ are used to calculate the timestep which is limited by variable Xi . The corresponding limited timei is calculated as follows: step Δtlim

2. Overview of the ATHROC code 2.1. Code structure and modelling method

i Δtlim =

The ATHROC code employs the idea of modular code development and Fig. 1 depicts the code structure. The modular code structure contributes to the convenience of code maintenance and improvement. The code includes two input files. One is the data file for setting up containment model. The other is accident-specific file for inputs like setting intervention or changing sensitive parameter. Changes to any data file inputs can be made in the accident-specific file. Therefore, the input flexibility helps to reduce the time to set up containment model and reduce input errors. Moreover, the structure of data file is not based on normally used numbered input style. Most of the input variables in the data file are given with a variable name followed by a required input value. All the variables are named according to a naming rule that is easy to be recognized. The detailed meanings of all the variable names are given in the user’s manual. Therefore, the input file is more readable and the code users can easily set up the containment model in the data file. The code treats the containment as a network of interconnected cells, in which the lumped parameter method is used. Each cell can represent one or more actual compartments in the containment, although sometimes the user may subdivide an actual compartment into several virtual cells to model some refined phenomena within the compartment, such as gas stratification. The cells communicate with each other through junction flow and/or conduction heat transfer in the heat structures. It is flexible to nodalize the containment and specify junctions, so many containment types can be modeled, especially the advanced light water reactor (ALWR) containment. The input flexibility also makes it easy to modelling some experimental facilities for code assessments. ATHROC uses two types of bulk fluids to model the transient behaviors of multiphase flow in the containment: the atmosphere bulk fluid, consisting of steam, non-condensable gases and homogeneously dispersed liquid droplets in the atmosphere, and the pool bulk fluid. Lumped parameter method, which homogenizes the properties of a

Xi ·FXi Xi̇

(1)

The final limited time-step Δtlim is set as the minimum value of all the limited time-steps. Δtlim was usually different from Δt0 , so Δt0 in the calculation of rates-of-change is adjusted to Δt ′0 and iteration calculations are performed until new Δt ′lim is within acceptable range of Δt ′0 . Finally the minimum value of Δt ′lim and Δt ′0 is set as the integration time-step. Integration calculation is performed to advance all the thermal-hydraulic variables to the next time step. This method adjusts the time step according to the intensity of transient process, thus, the computation speed is significantly accelerated. 3. General description of the code models 3.1. Thermodynamics model in one region The fluid properties within a given control volume were homogenized in the code. However, water and gas phases in a cell were considered to be thermal non-equilibrium, which allowed them to have

Fig. 1. Structure of ATHROC code. 77

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Since the exact form of velocity profile in the radius direction was not important in determining the physic of rising plume, uniform velocity and density profile was assumed at each height z. Then the control equations that describe the mass conservation of released light gas, mixed plume gases and momentum conservation of mixed plume gases can be written as

Start tart Read and process input files Initialize all the thermalhydraulic variables Initial run and reserve reserve Calculate the rate-of-change rate-of-change of variables by assuming all the comtributions

πρl uR2 = ṁ l

(6)

d (πρuR2) = 2πRue ρh dz

(7)

d (πρu2R2) = π (ρh −ρ) gR2 dz

(8)

where ṁ l is the mass rate of released light gas. The density of mixed gas in the plume, ρ , can be written as the sum of heavy entraining gas and released light gas,

Chose an appropriate time step and integrate the thermal hydraulic variables

ρ = ρh + ρl

(9)

The entraining velocity ue can be calculated by the following wellknown entrainment correlation (Ricou and Spalding, 1961):

Write different different output files

1

Last time step?

ρ 2 ue = E0 ⎛⎜ ⎞⎟ u ρ ⎝ ∞⎠

No

(10)

where E0 was taken as a constant value, 0.1. According to the state equation of ideal gas,

Yes

Stop Fig. 2. The flow chart of code computation.

P∞ = ρ∞

R R R T = ρl T + ρh T Mh Ml Mh

(11)

following equation can be obtained,

different temperatures. This assumption was appropriate for large cells containing little water, such as a containment compartment. Gas temperature of each compartment was computed from the total internal energy and the mass of steam and non-condensable gases:

ρ ρ ρ∞ = l + h Mh Ml Mh

(12)

The initial conditions for above equations are

Ug = ms us (Tg ,vs ) + Tg

∑

mi c vi

i

u = u 0,

(2)

∑ i

R = R0

at

z=0

ni RTg Vg

(3)

where Ps was obtained from steam table given gas temperature Tg and density of steam ρs . When the water pool was superheated due to heat addition or depressurization, the steam generation rate from the pool surface was determined from the degree of pool superheat, as follows:

Ws =

m w (h w−h w,sat ) (hs,sat −h w,sat )Δt

(4)

In contrast, when the steam was supersaturated, the rainout rate of the steam was determined from the degree of steam saturation, as follows:

Wrain =

ms (us,sat −us ) (us,sat −u w,sat )Δt

(13)

Finally, the dependent variables u,ρ,R can be obtained as functions of the plume height, z, then the flow rate of the rising plume is obtained

The total gas pressure was calculated by Gibbs-Dalton law assuming that the non-condensables were ideal gases:

P = Ps + Pnc = Ps +

ρ = ρl,0 = ρ0 ,

(5)

3.2. Buoyant plume model During an accident, light gaseous effluent such as hydrogen may be released from the primary system into the containment and rises like a plume, which would be accelerate by buoyancy and slowed down by the entrained surrounding gas. When the plume reaches the celling junction, mass exchange of gas between the two adjacent compartments occurs at two opposite flow directions. The complex process of this physical phenomenon is shown in Fig. 3.

Fig. 3. Sketch of buoyancy plume movement. 78

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

by WP = πR2ρu . The plume exiting a junction is restricted because of the finite opening area. Half of the junction opening area is assumed for the exiting plume and the other half for the incoming exchange flow,

Rout =

⎛ −Kij |Wij |Wij + ⎜Piex + ⎝ ⎛ − P jex + ⎜ ⎝

Ajunction (14)

2π

2

2

Ik (i)

k∈ i

Ik (i)

k∈ j

∂Pj ∂Wk

∂Pi P ⎞ − ρi gZ Wk + i δTi ⎟ e Pi ij,i ∂Wk Ti ⎠

Wk +

ρ

⎞ − j gZ δTj e Pj ij,j = 0 ⎟ Tj ⎠ Pj

(20)

For a containment system with J junctions and N cells, Eq. (20) represents only the J flow equations. However, there are J unknown flows and N unknown temperature changes in the equation set. Another N temperature change equations including δT should be introduced as described below. The change of atmosphere energy in cell i over the time step can be expressed in terms of “extrapolated” energy change and energy change due to atmosphere flow and cell temperature change over the current time step,

In order to obtain the flow rate of plume that leaves a circular 1 opening, a normal velocity with standard deviation σ = 2 R is assumed. Then the flow rate is calculated by

Wout = WP (1−e−2Rout / R )

∑

∑

(15)

3.3. Flow model between connected cells

ex

∂U δUi = ⎛ i ⎞ Δt + ⎝ ∂t ⎠

In this work, the flow between connected cells was modeled by junctions, which can transport either the atmosphere bulk fluid or the pool bulk fluid or both. The flow of atmosphere and pool was solved by independent momentum equations that neglects the shear and void fraction effects within the junctions. However, effects related to the coverage of the inlet or outlet of a junction by water pool were taken into account, such as gas-pool equilibration and liquid head effects, blockage of junction inlets by the pool, liquid head terms, and scrubbing effects.

∑

Ik (i) Wk (hj + c pj δTj )Δt + HAi δTi

k∈ i

(21)

∂U where ( ∂t i )ex

is the rate of change in Ui excluding the effect of flows. HAi is product of heat transfer coefficient and area, summed over all heat sinks in cell i. Alternatively,

δUi = mi δui + (ui + δui ) δmi = mi c vi δTi ex

⎡ ∂m + (ui + c vi δTi ) ⎢ ⎛ i ⎞ Δt + ∂t ⎠ ⎣⎝

∑

⎤ Ik (i) Wk Δt⎥ ⎦

3.3.1. Gas flow model For atmosphere bulk fluid, the flow was mainly driven by the pressure difference between cells. However, the cell pressure depended not only on the atmosphere flow, but also on the variations of the atmosphere temperature, which was associated with the heat and mass transfer to heat sinks. Consequently, an implicit numerical method was used where flow rates were calculated based on end-of-time step cell pressures and atmosphere temperatures. Solving flow rates based on current pressures would be numerically stiff and required very small time steps. The quasi-steady flow model was used assuming that the fluid inertia was negligible. The momentum conservation equations are referenced to Murata et al. (1997) as follows,

Equations (20) and (23) form a J + N non-linear system of equations with J + N unknown variables, which can be written as the following matrix form and solved by Newton-Raphson method.

K |Wij |Wij = ΔP

F (x ) = 0

ex

⎛ ∂Ui ⎞ Δt + ⎝ ∂t ⎠

(18)

The implicit pressure Pi or Pj can be expressed in terms of extrapolated pressure and pressure changes due to flow and atmosphere temperature change:

∑ k ∈ i

Ik (i)

∂Pi P Wk + i δTi ∂Wk Ti

∑ k∈ i

⎤ Ik (i) Wk Δt⎥ ⎦

(23)

(24)

3.3.2. Countercurrent flow model Combined buoyancy-driven exchange flow and forced flow through an opening at a ceiling or floor was considered when the two ventconnected cell were filled with fluids of different density together with an unstable configuration, e.g., density at the top was larger than the bottom (Cooper, 1995). In case of small pressure difference, bidirectional exchange flow may occur in the opening. But for larger pressure difference, the content of buoyancy-driven exchange flow would be prevented and only unidirectional forced flow occurs. A simple case of such combined flow through an opening at a vertical wall is shown in Fig. 4. Fluid is supplied to cell 1 and extracted from cell 2 at a volume flow rate Q0 . If Q0 increase from zero to a large value, the flow configuration through the opening goes through pure buoyancy-driven countercurrent flow, combined flow and pure forced flow. QF −QBF is equal to Q0 for all cases, and the relation between flow configuration and Q0,QF ,QBF is shown in Table 1. “Flooding flow” is the critical flow regime between combined flow and pure forced flow, and its corresponding back forward flow equals to zero. Epstein and Kenton (1989) performed an experiment using water and brine as a substitute for gases with different densities. Pure countercurrent flow rate, flooding flow rate and back flow rate were obtained for vertical walls and horizontal ceilings. The correlations were used by the code (Epstein and Kenton, 1989).

then following general form could be derived for the flow through a junction with the central height of Zij,i and Zij,j at the two ends,

Pi = Piex +

Ik (i) Wk (hj + c pj δTj )Δt + HAi δTi ex

(17)

ρj ρ −( ) gZ −( i ) gZij,i Pi −Pj e Pj ij,j

∑ k∈ i

⎡ ∂m = mi c vi δTi + (ui + c vi δTi ) ⎢ ⎛ i ⎞ Δt + ∂t ⎠ ⎣⎝

where K is the junction flow resistance determined by flow regime. Since the pressure in a control volume can be described by

K |Wij |Wij = Pi e

(22)

Combining Eqs. (21) and (22) yields the following temperature change equations:

(16)

dP = −ρ (P ) g dZ

k∈ i

(19)

where the “extrapolated” pressure, Piex, is defined as the end-of-time step pressure if no inter-cell flow exits and if the compartment gas temperature is constant. δTi is the temperature change of cell i due to atmosphere flow and heat transfer to heat sinks over the current time step. Ik (i) is the identification number of junction flow Wk , which has a value of 1 if cell i is the receiver cell of the flow and a value of -1 inversely. Substituting Eq. (19) into Eq. (18) gives 79

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

gaseous bulk fluid and a slab. Mechanistic heat and mass transfer analogy (HMTA) methodology were developed to treat this process. In addition, optional semi-empirical correlations such as Tagami and

Pj

Pi

Zj

Zi

ZWi

ZWj

Fig. 4. Sketch of countercurrent flow.

(a) One side completely uncovered, the other partially covered

Table 1 Flow configuration identification.

Pj

Flow configuration

Values of Q0,QF ,QBF

Pure buoyancy-driven exchange flow Combined flow

Q0 = 0,QF = QBF > 0 0 < Q0 < Qflood , QF > QBF > 0

Flooding flow

Q0 = Qflood , QF > QBF = 0

Pure forced flow

Q0 > Qflood , QF > QBF = 0

Pi

ZWj

Zj

3.3.3. Water flow model The treatment for water flow between connected cells was very similar to the atmosphere flow. However, the junction ends could be uncovered in the water flow when the flow calculations were converged. As in the atmosphere flow model, the water flow could be solved by a set of momentum equations describing the flow in each junction. For a normal junction, the flow equation was written as

|Wij |Wij = 2Aij,0 Aij,n ρi CD2 ΔP

ZWi

(b) Both sides partially covered Pj

Pi

(25) ZWi

ΔP can be expressed in terms of water gravitational head and gas pressure difference, that is

ΔP = g {ρwi (Z wi−Zi )−ρwj (Z wj−Zj )} + ΔPL + (Pi−Pj )

Zi

Zi Zj

ZWj

(26)

(c) One side completely uncovered, the other completely covered

where ΔPL is the water gravitational head due to elevation difference of the two junction ends. For horizontal junctions, ΔPL is equal to zero. The five horizontal flow regimes considered in the code are shown in Fig. 5. Note that cell i and cell j were just pre-assigned donor and receiver cells, and they were replaced with others when a negative flow was calculated for the junction. When at least one side of the junction was not completely covered, gas flow through the uncovered part of the junction is assumed sufficient to balance Pi and Pj. In this case, Pi - Pj can be removed from Eq. (26).

Pj

Pi

Zj

ZWj

ZWi

Zi

3.4. Heat transfer model for passive heat sinks

(d) One side partially covered, the other completely covered ATHROC treated a heat sink as a slab, assigned properties to each slab, nodalized the slab, and specified boundary heat fluxes that are imposed to the surfaces of the slab. The heat transfer coefficients, steam condensation rate, and energy exchange rate for the heat sinks were calculated. Thermal resistances such as paint layers, metal liners, gas gaps were considered for conduction heat transfer. The heat transfer mode consisted of radiation, natural convection, forced and mixed convection. In case of a primary system blowdown into the containment, vapor condensation process may occur on the heat sinks. ATHROC could simulate the condensation film flow on the surfaces of contiguous walls which extends from one control volume to another volume. The contiguous walls were subdivided into several connected heat sinks. Heat and mass transfer were coupled at the interface between the

Pj

Pi

ZWi ZWj

Zj

Zi

(e) Both sides completely covered Fig. 5. horizontal flow regimes. 80

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Uchida correlations were also available for the design basis analysis. The radiation heat transfer coefficient was given by

hrad =

σ (Tg4−Tw4 )

(

1 Eg

+

) (T −T )

1 −1 Ew

g

w

(27)

Note that when the temperature difference between gas and wall surface, Tg – Tw, was less than or equal to 1 K, 'Hopital's rule was applied and the resulting heat transfer coefficient was

hrad =

4σTg3

(

1 Eg

+

)

1 −1 Ew

(28)

The natural convective heat transfer coefficient hnc was given by

Nu = hnc L/ λ g = a (Ra)n

Fig. 6. sketch of energy balance of the film.

(29) For a contiguous wall which was subdivided by several connected heat sinks, film will flow along the wall that extended from one control volume to another. In this case, the condensation heat transfer coefficient was given by Seban’ correlation (Seban and Chun, 1971):

where a and n are constants associated with surface configuration and flow regime. Ra, Gr, Pr can be calculated as follows,

Ra = Gr Pr, Gr =

g|ρg −ρf |L3 Tg νg2

, Pr =

μg cpg λg

(30)

hcd = 0.606 ⎛⎜ ⎝

The forced convective heat transfer coefficient was given by Kreith et al. (2012) as follows,

⎧ Nu = 2; 1/3 ⎪ ⎪ Nu = 1.86(RePr)1/3 DH ; L ⎨ ln (Nu) = a ln(Re) + b ; 1 1 ⎪ ⎪ Nu = 0.23Re0.8Pr 0.3; ⎩

hm =

hc Pr 0.66 ⎛ ⎞ c pg Pam Mw ⎝ Sc ⎠

q′atm = (hc + hrad )(Tg−Tf ) + hm Mw [Ps−Psat (Tf )] hsat (Ps )

qhs = qatm; still film qhs + qout = q′atm + qin;

where L is the dimensionless length defined by 1/3

2 ⎛ ρl g sinθ ⎞ ⎟ Pr[hfg + 0.68Cl (Tf −Tw )] ⎝ μl2 ⎠ ⎜

1 ρg

1− ρ

l

film

(41)

During a LOCA or MSLB accident, spray system is essential to prevent the containment from over-pressurization. The spray heat removal model computed the mass and energy exchange of atmosphere and spray droplets by means of evaporation, condensation and convective heat transfer. Typically, the droplets entered at a cold temperature below the dew point of steam in the atmosphere. Steam was condensed on the droplets, which were heated by convective and latent heat transfer. The droplets began to evaporate when its temperature exceeded the dew point. They reached a wet bulb temperature when the convective heat transfer was balanced with evaporative cooling. The evaporating droplet continued to drop down until it was entirely evaporated or it reached the water pool surface in the cell. The heat and mass transfer process for a droplet could be expressed by the following conservation equations (Kreith et al., 2012):

(34)

4L (Tf −Tw ) Cl

flowing

3.5. Spray heat removal model

1

L=

(40)

Finally, the temperature of film surface, Tf, was iteratively determined by following energy balance equations,

The heat transfer through a condensation film and into the heat sink was determined by the orientation of wall surface. For a vertical or inclined wall, the heat transfer coefficient, hcd, for both the turbulent and the laminar film condensation was calculated based on the work of Rohsenow et al. (1998):

hcd

(39)

The mass conservation equation for the flowing film is

Wout = Win + Δmcd = Win + hm AMw [Ps−Psat (Tf )]

(33)

ρ 2 g sinθ 3 = λl ⎜⎛ l 2 ⎟⎞ 10 f (ln(L)) ⎝ μl ⎠

(38)

For flowing film, the heat flux from the gas to the film (not to the heat sink surface) could be calculated as follows:

|P−Ps−(P−Psat (Tf ))| P − Ps | sat (Tf )

(37)

qatm = (hc + hrad )(Tg−Tf ) + hm Mw [Ps−Psat (Tf )](hs−hcon )

(32)

ln| P − P

(36)

Alternatively, for still film, the heat flux from atmosphere to heat sink could be calculated by

where Pam is the log mean noncondensable gas pressure in the free gas phase at the condensing surface,

Pam =

Re−l 0.22

⎠

qhs = hcd (Tf −Tw )

(31)

Note that the constants a1 and b1 in the equation for transition zone from laminar to turbulent flows, 2000 < Re ⩽ 6000,was interpolated from laminar and turbulent correlations. The heat transfer coefficients provided above was used to evaluate sensible heat transfer to the wall. However, when condensation or evaporation occurs, additional latent heat transfer should be taken into accounted. In this case, heat and mass transfer process were coupled at the condensation film surface. Through the interface energy balance equation, the interface temperature Tf which drives the interface processes could be determined. The energy balance equations were discussed as follows. The mass transfer coefficient was evaluated using Reynold's analogy between heat and mass transfer (Collier and Thome, 1994; Eckert and Drake Jr, 1987):

Sc 0.66 Sh = Nu ⎛ ⎞ ; ⎝ Pr ⎠

0.33

⎟

Wl is the average flow rate of the film, Pl is the userwhere input flow perimeter for each heat sink. Fig. 6 shows the sketch of energy balance of the still film and flowing film. The heat flux through both films could be calculated as

10 < Re⩽2000 2000 < Re⩽6000 Re > 6000

μl2

4W Re= P μl , l l

Re⩽10

( )

λl3 ρl2 gsinθ ⎞

dmd = hm A [Pst −Psat (Td )] dt

(35) 81

(42)

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

md c p

dTd = hc A (Tg−Td ) + hm Ahfg (Pst −Psat (Td )) dt

(43)

Nu = 2.0 + 0.6Re1/2Pr1/3, Sh = 2.0 + 0.6Re1/2Sc1/3

Wet part

170 cm

where the convective heat transfer coefficient hc and mass transfer coefficient hm could be determined from the following correlations that used Reynold analogy method (Ranz and Marshall, 1952): (44)

Tg RPam d Mw Ddiff P

, Sc =

μg ρg Ddiff

, Re =

ρg Ud d μg

80 cm

Sh = hm

(45)

Dry part

439.2 cm

The dimensionless numbers in Eqs. (44) were defined as

The droplet velocity Ud was obtained by balancing the buoyancy force with the drag force, that is 1/2

⎡ 4gd (ρd −ρg ) ⎤ Ud = ⎢ 3ρg CD ⎥ ⎣ ⎦

(46) 180 cm

4. Code assessments A good code validation matrix consists of all the important phenomena occurring in the containment. In order to set up the code validation matrix for ATHROC, two NEA/CSNI reports (Chin et al., 2014; Karwat et al., 1999) were intensively reviewed, which discussed PIRT (Phenomena Identification and Ranking Table) on the containment thermal hydraulic phenomena and defined a basic set of available experiments for code validation. In this paper, typical separated effects and integral effects tests like Phebus FPT0 test, JAERI tests, CVTR tests, and NUPEC M-7-1 test were simulated to assess the code reliability in predicting containment thermal hydraulics. Table 2 shows some important containment phenomena addressed in these tests. Details can be referenced to the above two NEA/CSNI reports.

Sump Fig. 7. Phebus containment vessel.

temperature controlled surfaces. The injected steam could only condense on the wet part. The outer containment walls were slightly overheated to prevent any steam condensation. The lower containment vessel structure included an elliptical bottom and a 0.1 m3 sump with a diameter of 0.58 m in order to model atmosphere water-surface. Saturated steam was injected into the lower region of containment vessel at injection rates that varied from 0.0005 to 0.003 kg/s, as shown in Fig. 8. The steam enthalpy was 2746.76 J/kg. This test was a typical separate effects test for natural convection condensation. It was simulated with the ATHROC code using a single cell to represent the vessel atmosphere. Sensible heat transfer by convection and radiation was calculated in the code. Natural convective heat transfer correlation for vertical walls was used to predict the Nusselt number at the condenser surface. Condensation in the code was modeled using the heat and mass transfer analogy method. Heat transfer to the containment vessel was only by sensible heat transfer (convection and radiation) which was significantly less than the latent heat transfer that occurs on the condensers. Because of the relatively high wall temperatures and steam injection temperature, the bulk atmosphere was superheated throughout the test. Fig. 9 shows the comparison of measured and calculated pressure during the transient period. A similar trend was predicted, while the

4.1. Phebus FPT0 test The PIRT information provided by Karwat et al. (1999) indicated that heat and mass heat transfer to the heat sinks were highly ranked phenomena. So the benchmark of Phebus FPT0 test was selected to assess the condensation model. The Phebus FP (Fission Products) international research program was conducted between 1988 and 2010 and the FPT0 test was realized in 1993. In this study, only the thermal hydraulic behavior of the steam/gas mixture in the containment vessel in the FPT0 test was considered. As shown in Fig. 7, the stainless steel containment vessel had a volume of 10 m3 and the condensation heat transfer on the heat sinks were studied by installing three vertical condensers centrally located in the upper region of the containment. Each condenser was subdivided into a 170 cm long wet part and an 80 cm long dry part with Table 2 Important containment phenomena addressed in these tests. Component Atmosphere Pressurization Mixing

Transport Structure interior Heat transfer Structure surface Mass transfer (condensation/evaporation)

Phebus FPT0

JAERI

CVTR

NUPEC M-7-1

√

√ √

√ √ √ √ √ √ √

√ √ √ √ √ √ √

√

√

√ √

√

Phenomena

Gas compression/expansion Spray mass and heat exchange Jet plume gas interaction/entrainment Buoyancy/stratification Spay dynamics Buoyancy Form and friction losses

√ √

1-D heat conduction √

Free convection Forced convection

82

√

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

water is injected, the vessel contains a saturated air/steam mixture at 3.5 bar. The vessel steel walls are insulated and pre-heated to prevent steam condensation on the surfaces. In our simulation, single cell model was used for the two tests. For PHS-6, no spray water was assumed to contact the vertical walls of the vessel. For PHS-1, some injected spray water was assumed to contact the vertical wall and therefore vaporize, thus reducing the spray cooling rate. A default spray droplet diameter of 0.001 m was used. The vertical walls were modeled with 9 continuous heat sinks to track film flow on the surfaces. Fig. 13 shows the comparison of predicted pressure and experimental data for PHS-6 test. The good agreement confirmed the spray heat and mass transfer model in the code. The comparison of calculated gas temperature and data is shown in Fig. 14, where the temperature data represent the uniform temperature at different elevations. A slight increase in calculated temperature at 6000 s was observed, indicating that the thermal condition transfers from a saturated state to a slightly superheated state. However, the superheated condition sustained a short time and quickly transfers to a saturated condition. It should be noted that superheating was not a significant concern for most containment loads analysis, because sensible heat transfer due to superheating was smaller than latent heat transfer associated with condensation. Small superheating only resulted in minor pressure variations. Fig. 15 shows the comparison of predicted pressure and data for PHS-1 test, in which 10 percent of injected spray water was assumed to contact the vertical walls and resulted in film flow on the surfaces. The film flow was modeled with film tracking model in the code. The predicted rate of pressure decrease was slightly larger than the data, but the relative error was less than 10%. Referring to the good agreement of pressure in PHS-6 test, the small difference may be caused by the inaccuracy of film tracking model. Comparisons of calculated temperature and measured data at different elevations are shown Fig. 16. Predicted temperature agreed well with the data at locations below the spray injection. However, large deviations were observed for the region near the top of the vessel where the assumption of uniform mixing was invalid.

Fig. 8. Steam and hydrogen injection rates during FPT0 test.

Fig. 9. Comparison of pressure between data and AHROC results.

4.3. CVTR tests calculated pressure was slightly higher than the experimental data. The maximum deviation from data was less than 5%. Pressure response was just an indirect indicator to validate the condensation model, because it depended not only on the total energy of the gas, which was affected by the condensation rate, but also on the thermal-dynamic state of gas mixture. Condensation rate was a direct indicator for the accuracy of condensation model, and Fig. 10 shows the comparison of calculated condensation rate with the experimental data. The very good agreement indicated that the natural convective condensation model in the code was accurate for FPT0 test conditions. Fig. 11 shows the comparison of the measured and calculated saturation and superheated temperatures. The experimental data of saturation temperature were determined from total pressure data and relative humidity (1 for FPT0 test). As shown in the results, calculated saturation temperature and superheated temperature both agreed well with the experimental data, indicating that the heat and mass transfer model in ATHROC was applicable for superheated conditions.

The CVTR facility is a decommissioned reactor containment building. A series of DBA simulation tests were conducted in CVTR in the late 1960’s to provide experimental data for code validation. The test facility and the series of experiments performed in this facility can be found in Schmitt et al. (1970). Fig. 17 shows the test configuration

4.2. JAERI tests JAERI tests are essentially separate effects tests for sprays cooling. They were conducted in a 700 m3 steel vessel (20 m high, 7 m in diameter) during the late 1970s in Japan. Two tests, PHS- 1 and PHS-6, were considered in this paper. Fig. 12 is a sketch of the configuration for the two tests. PHS-6 is a single nozzle test with an 18 m spray height, while PHS-1 is a six nozzles test with a 15 m spray height. Before spray

Fig. 10. Comparison of condensation rates between data and AHROC results. 83

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Fig. 11. Comparison of saturation and superheated gas temperature between data and AHROC results.

Fig. 14. Comparison of predicted and measured gas temperature for PHS-6 test.

20 m PHS-6

18 m

PHS-1

15 m

Fig. 15. Comparison of predicted and measured pressure for PHS-1 test.

sump Fig. 12. JAERI vessel and test configuration for PHS-6 and PHS-1.

Fig. 16. Comparison of predicted and measured gas temperature for PHS-1 test.

and code nodalization for CVTR facility. The total free volume of the containment is about 6426 m3, and the junction between operational region and intermediate region is 30.6 m2. Spray nozzles are uniformly spaced at the containment bend line and the average spray droplet diameter is suggested to be 1 mm. The steam was injected at an

Fig. 13. Comparison of predicted and measured pressure for PHS-6 test.

84

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Fig. 17. Test configuration and code nodalization.

elevation slightly above the operation deck. Accordingly, the geometry of the containment structure was appropriately modeled by 19 cells. The heat sink input was based on “best-estimate” concrete areas given in the report (Schmitt et al., 1970). Three simulated design basis tests, test 3, 4 and 5 were conducted, which replicates the postulated design basis “intermediate” size MSLB events. For all these tests, the steam release data were nearly identical, as indicated in Table 3. Sprays were not activated for Test 3, but were activated at about 200 s with a flow rate of 0.0183 m3/s and 0.0316 m3/ s for Test 4 and Test 5, respectively. Fig. 18 shows the comparison of calculated dome gas pressure and experimental data. The measured peak pressures during the blowdown period for the three tests were 204 kPa, 229 kPa, 228 kPa, respectively, which were slightly larger than the predicted values. The reason for this diffenrece may be that in the initial pressurization process, momentum induced mixing and forced convection heat and mass transfer play an important role in the containment response, but these phenomena were not well considered in the code. From a long term respective (1000 s3000 s), pressures decreased at a rate approximately equals to the experimental data for all the three tests. Fig. 19 shows the short-term gas temperature in the dome. The peak temperature and the decrease of temperature during cooldown were well predicted by the code. Underestimated peak pressure and well predicted peak temperature indicated that the calculated containment condition at the end of blowdown period has a more degree of superheat than experimental data. The decrease in gas pressure during the cooldown period of case 4 and case 5 (saturated condition) were well predicted, but the predicted decrease

Fig. 18. Comparison of predicted and measured dome pressure.

rate of gas pressure for case 3 (superheated condition) was smaller than measured. These indicated that the heat and mass transfer model in the code was more applicable to saturated conditions. Fig. 20 shows the long-term gas temperature in the dome. The containment conditions were driven by natural convection process in the long term. The deviation between calculated and measured temperature was small, which provided a strong evidence that the natural convection heat and

Table 3 Steam Data. Test No.

Steam injection time of duration (s)

Superheated steam flow rate (kg/s)

Saturated water flow rate (kg/s)

Saturated steam flow rate (kg/s)

Total flow (kg)

Average enthalpy (J/kg)

3 4 5

166.4 174.8 173.1

40.32 39.06 40.32

7.56 8.19 5.04

47.88 47.25 45.36

7495.6 7865.7 7428.0

2779528.18 2781854.14 2784180.10

85

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

chimney loops (Comp. 20, 21, 23, 24). The spray water in each lower compartment was supplied with a proportion of spray from the dome. The spray was assumed to concentrate on the center of the containment so that no spray water would impinge on the containment wall. Any spray that did not fall into the lower compartments was collected in each compartment floor and eventually drained into the water storage tank (Comp. 26). The spray was assumed to be uniformed over the cross section of each compartment. The steam/helium mixture was injected into the SG foundation compartment (Comp. 8). The comparison of predicted and measured pressure in the dome is shown in Fig. 22. The general trend was well predicted. The pressure history could be divided into 3 phases before 40 min. Phase 1 and Phase 2 were the 30-min period when sprays were operating. During Phase 1, pressure dropped due to spray water cooling on the atmosphere. The rate at which the pressure dropped was well predicted by ATHROC. During Phase 2, pressure increased as the hot gas mixture was injected and accumulated in the containment. During Phase 3, pressure increased at a larger rate due to the termination of water spray. The rate at which pressure increased was over-predicted during phase 2 and phase 3. The reason was that the code failed to distribute some of the spray water to the internal heat sinks, such as SG room walls, resulting in insufficient cooling of the heat sinks and large quantities of returning heat from heat sinks to surrounding atmosphere. On an absolute basis, the predicted pressure at 30 min and 40 min was approximately 2.3% and 3.8% greater than the data respectively. Fig. 23 shows the temperatures in the SG foundation compartment (Compt. 8) and all of the rooms directly above it (Compt. 15, 21, 25). Fig. 24 shows the temperature in another vertical column of SG rooms. Fig. 25 shows the temperatures in a vertical column of outer rooms: the lower general compartment (Compt. 4) and the rooms directly above it. Fig. 26 shows the temperatures in a number of central rooms, including In Core Chase room (Compt. 1), C/V sump pump room (Compt. 2) and cavity (Compt. 19). Fig. 27 shows the results for pressurizer rooms, including a dead-end compartment (Compt. 16). In general, the trends of the predicted temperatures agreed well with the experimental data. Moreover, the temperature difference between lower compartments and the dome was also well predicted. The maximum deviation from measured gas temperatures was less than 10% on an absolute basis. The temperature history was characterized by a continuous decrease during spray period, though there was, at the beginning, a significant temperature increase in the gas mixture injection compartment (Compt. 8) and two compartments directly above it (Compt. 15 and 21). Relatively large deviation between predicted temperature and experimental data was presented in the Comp. 1, which has a relatively closed, complex geometry and low elevation inside the containment. There were also some differences between predicted and measured temperature in the pressurizer compartments (Comp. 16 and Comp. 22). The deviation may be caused by the lack of wall cooling by spray water distributed to the outer surface of Comp. 16 and 22. Figs. 28–30 show the results of predicted and measured helium

Fig. 19. Short-term gas temperature variation in the dome.

Fig. 20. Long-term gas temperature variation in the dome.

mass transfer model in the code was reliable. 4.4. NUPEC M-7-1 test NUPEC M-7-1 test, which was selected as an international standard problem (ISP-35), was simulated using ATHROC code to assess the code's capability for beyond design basis accident analysis where hydrogen gas was released into the containment. The effects of mixing from natural convection as a result of density differences, forced convection as a result of the steam release and water sprays were tested. The test conditions are shown in Table 4. The NUPEC facility is a 1/4 linearly scaled PWR dry, insulated steel containment vessel with a free volume of 1300 m3. There are 25 subcompartments separated by steel partitions. The large dome occupied 70% of the total free volume. The containment has 3 main floors located at elevations of 3200, 5425 and 7325 mm, respectively. 21 water spray nozzles are installed in the top region of the dome. Inner and containment walls are made of carbon steel. In addition, the containment wall and the first floor are covered with insulation. Fig. 21 shows the 26-compartment nodalization scheme of the NUPEC facility for ATHROC code. A total of 67 junctions connecting these compartments were established. In addition, inner walls, containment walls, floors, insulators were modeled with 62 heat sinks. Water sprays were injected into the dome (Comp. 25) and lower SG

Table 4 NUPEC M-7-1 test conditions.

86

Item

Parameter

Initial pressure Initial temperature He temperature Steam temperature Steam flow rate He flow rate Gas injection duration Spray Water temperature Spray rate Spray duration

∼140 kPa ∼70 °C 15 °C 160 °C 0–0.03 kg/s 0.08–0.33 kg/s 30 min 40 °C 19.4 kg/s 30 min

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Fig. 24. Comparison of temperatures in Comp. 7, 14, 20, 25.

Fig. 21. Containment nodalization scheme.

Fig. 25. Comparison of temperatures in Comp. 3, 4, 12, 25.

Fig. 22. Atmosphere pressure in the dome.

Fig. 26. Comparison of temperatures in Comp. 1, 2, 19, 25.

concentrations in different compartments. The helium was fairly well mixed throughout the containment except for the vertical column of injection compartments (Comp. 8, 15, 21). This was primarily resulted from the mixing effect of water sprays, but was also due to the low elevation of gas injection and the relatively open geometry of the facility. The injection of hot gas mixture in Comp. 8 along with the rapid pressure drop in the dome resulted in a high gas flow to the dome.

Fig. 23. Comparison of temperatures in Comp. 8, 15, 21, 25. 87

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Fig. 27. Comparison of temperatures in Comp. 16, 22, 25.

Fig. 30. Comparison of helium concentrations in Comp. 16, 22, 25.

inaccuracy flow loss coefficients used in our simulation. Helium concentration in the other compartments was relatively well predicted by ATHROC. Even the pressurizer compartments, Comp. 16, 22, had a maximum deviation less than 16%. The helium source of Comp. 16 and 22 is from the dome, so these three compartments had the helium concentrations directly proportional to the height of the compartments. It should be noted that Comp. 16 was a dead-end compartment which was not easy to modelling precisely using lumped parameter code. However, ATHROC adopted counter-current flow model to deal with this problem, and the results for Comp. 16 and 22 illustrated the accuracy of this model. 5. Conclusions In this study, the home-made code ATHROC was developed to predict the long term thermal hydraulics behaviors (including hydrogen distribution) of China advanced PWR containment during design basis accidents or some core degrade accidents. General containment model and some supplemented models such as plume model and countercurrent flow model for sub-node physics were developed. Moreover, specially designed method was applied to speed up the calculation. ATHROC code assessment was carried out by simulating several separated effects tests and integral effects tests. Most of the results showed that the code could provide reasonable predictions. Some important findings and conclusions could be summarized as follows:

Fig. 28. Comparison of helium concentrations in Comp. 8, 15, 21, 25.

• Results for Phebus FPT0 Test showed that the natural convective •

Fig. 29. Comparison of helium concentrations in Comp. 7, 14, 20, 25.

•

Consequently, helium concentrations in Comp. 15 and 21 were close to that in Comp. 8. Before 15 min, linear increase in the helium concentration of Comp. 8 was consistent with the linear increase in the mass of injected helium. The predicted rate of increase and the peak value were both smaller than the data, but the time of peak value was well predicted. The maximum deviation from measured data was about 25% for Comp. 8, 15, 21.The large discrepancy may be caused by the 88

condensation model in the code could be applied to superheated conditions. However, the results for CVTR Tests indicated that the natural convective condensation model was more applicable to saturated conditions than superheated conditions. Results for JAERI Tests confirmed the spray heat and mass transfer model in the code. The default spray droplet size of 0.001 m was recommended. In the case of multiple nozzles spray, it was reasonable to assume that some of the spray water will contact the containment wall and then flow downward. Small superheating only resulted in minor pressure variations, thus it was not a significant concern for the containment loads analysis. Results for NUPEC M-7-1 Test proved the code’s capability of predicting hydrogen distribution under beyond design basis accident conditions, which involves the effects of mixing from natural convection and forced convection. In addition, meaningful predictions of temperature stratification can be calculated. However, the lack of modelling for allocating a portion of the spray water to the internal heat structures causes some simulation distortion for relative closed geometry.

Nuclear Engineering and Design 334 (2018) 75–89

Y. Chen et al.

Acknowledgements

management and related computer codes. Organisation for Economic Co-Operation and Development. Lin, C., 1979. CONTEMPT-LT/028. Pressure-Temperature Response, EG and G Idaho Inc, Idaho Falls, ID (United States). Lin, C., Economos, C., Lehner, J., Maise, G., Ng, K., 1986. CONTEMPT4/MOD6. PWR & BWR Containment Analysis, Nuclear Regulatory Commission, Washington, DC (United States). Motoki, Y., Naritomi, M., Tanaka, M., Nishio, G., Hashimoto, K., Kitani, S., 1983. Heat removal tests for pressurized water reactor containment spray by large-scale facility. Nucl. Technol. 63, 316–329. Müller, K., Toth, B., 2007. Final interpretation report of the PHEBUS test FPT0 (Bundle Aspects)–DG-JRC-IE. IE Technical Report October 2007 (EUR 23222 EN). Murata, K., Williams, D., Griffith, R., Gido, R., Tadios, E., Davis, F., Martinez, G., Washington, K., Tills, J., 1997. Code manual for CONTAIN 2.0: a computer code for nuclear reactor containment analysis. Nuclear Regulatory Commission, Washington, DC (United States). Div. of Systems Technology; Sandia National Labs., Albuquerque, NM (United States); Tills (Jack) and Associates, Inc., Albuquerque, NM (United States). Ranz, W., Marshall, W., 1952. Evaporation from drops. Chem. Eng. Prog 48, 141–146. Ricou, F.P., Spalding, D., 1961. Measurements of entrainment by axisymmetrical turbulent jets. J. Fluid Mech. 11, 21–32. Rohsenow, W.M., Hartnett, J.P., Cho, Y.I., 1998. Handbook of heat transfer. McGraw-Hill, New York. Royl, P., Rochholz, H., Breitung, W., Travis, J.R., Necker, G., 2000. Analysis of steam and hydrogen distributions with PAR mitigation in NPP containments. Nucl. Eng. Des. 202, 231–248. Schmitt, R.C., Bingham, G., Norberg, J., 1970. Simulated design basis accident tests of the carolinas virginia tube reactor containment. Final Report. Idaho Nuclear Corp, Idaho Falls. Seban, R., Chun, K., 1971. Heat transfer to evaporating liquid films. J. Heat Transfer 93, 391–396. Sehgal, B.R., 2011. Nuclear safety in light water reactors: severe accident phenomenology. Academic Press. Stamps, D.W., 1995. CONTAIN assessment of the NUPEC mixing experiments. SANDIA REPORT, SAND94-2880 UC-505, Albuquerque. Tills, J., Notafrancesco, A., Longmire, P., 2008. An Assessment of MELCOR 1.86: Design Basis Accident Tests of the Carolinas Virginia Tube Reactor (CVTR) Containment (Including Selected Separate Effects Tests). SAND2008-1224, Sandia National Laboratories, Albuquerque, New Mexico. Van Dorsselaere, J., Chatelard, P., Cranga, M., Guillard, G., Trégourès, N., Bosland, L., Brillant, G., Girault, N., Bentaib, A., Reinke, N., 2010. Validation status of the ASTEC integral code for severe accident simulation. Nucl. Technol. 170, 397–415. Van Dorsselaere, J., Pignet, S., Seropian, C., Montanelli, T., Giordano, P., Jacq, F., Schwinges, B., 2005. Development and assessment of ASTEC code for severe accident simulation. Xing, J., Song, D., Wu, Y., 2016. HPR1000: Advanced Pressurized Water Reactor with Active and Passive Safety. Engineering 2, 79–87.

The authors appreciate the financial support from National Natural Science Foundation of China (No. 11575137) and National Natural Science Foundation of China (No. 11605131). References Allelein, H.-J., Arndt, S., Klein-Heßling, W., Schwarz, S., Spengler, C., Weber, G., 2008. COCOSYS: status of development and validation of the German containment code system. Nucl. Eng. Des. 238, 872–889. Chin, Y.-S., Mathew, P., Glowa, G., Dickson, R., Liang, Z., Leitch, B., Barber, D., Vasic, A., Bentaib, A., Journeau, C., 2014. Containment Code Validation Matrix. Organisation for Economic Co-Operation and Development. Clément, B., Hanniet-Girault, N., Repetto, G., Jacquemain, D., Jones, A.V., Kissane, M.P., von der Hardt, P., 2003. LWR severe accident simulation: synthesis of the results and interpretation of the first Phebus FP experiment FPT0. Nucl. Eng. Des. 226, 5–82. Collier, J.G., Thome, J.R., 1994. Convective boiling and condensation. Clarendon Press. Cooper, L.Y., 1995. Combined buoyancy and pressure-driven flow through a shallow, horizontal, circular vent. J. Heat Transfer 117, 659–667. Eckert, E.R.G., Drake Jr, R.M., 1987. Analysis of heat and mass transfer. Epstein, M., Kenton, M.A., 1989. Combined Natural Convection and Forced Flow Through Small Openings in a Horizontal Partition, With Special Reference to Flows in Multicompartment Enclosures. J. Heat Transfer 111, 980–987. Gauntt, R., Cole, R., Erickson, C., Gido, R., Gasser, R., Rodriguez, S., Young, M., 2005. MELCOR computer code manuals. Sandia National Laboratories, NUREG/CR, pp. 6119. George, T.L., Wiles, L., Claybrook, S., Wheeler, C., McElroy, J., 2001. GOTHIC Containment Analysis Package Technical Manual. Numerical Applications Inc. Gyenes, G., Ammirabile, L., 2011. Containment analysis on the PHEBUS FPT-0, FPT-1 and FPT-2 experiments. Nucl. Eng. Des. 241, 854–864. Hong, S.J., Choo, Y.J., Hwang, S.H., Lee, B.C., Ha, S.J., 2015. Development of CAP code for nuclear power plant containment: Lumped model. Nucl. Eng. Des. 291, 47–63. Karwat, H., Bardelay, J., Hashimoto, T., Koroll, G.W., Krause, M., L'Heriteau, J.-P., Lundstroem, P., Notafrancesco, A., Royl, P., Schwinges, B., 1999. SOAR on Containment Thermal-hydraulics and Hydrogen Distribution-Prepared by an OECD/ NEA Group of Experts. Organisation for Economic Co-Operation and DevelopmentNuclear Energy Agency. Klein-Heßling, W., Arndt, S., Weber, G., 2000. COCOSYS V 1.2 User Manual. GRS, July. Kreith, F., Manglik, R.M., Bohn, M.S., 2012. Principles of heat transfer. Cengage learning. Layly, V.D., Spitz, P., Tirini, S., Mailliat, A., 1996. Analysis of the Phebus FPT0 containment thermal hydraulics with the Jericho and Trio-VF codes. Nucl. Eng. Des. 166, 413–426. Liang, Z., Chan, C., Sonnenkalb, M., Bentaib, A., Malet, J., Sangiorgi, M., Gryffroy, D., Gyepi-Garbrah, S., Duspiva, J., Sevon, T., 2014. Status report on hydrogen

89