1 Introduction

LOCA (Loss-Of-Coolant Accident) in the primary loop and MSLB (Main Steam Line Break) accidents are separately considered as design-basis accidents in reactor safety analysis. According to analysis of the two accident modes, the most serious LOCA-guillotine main coolant pipelines rupture for VVER1000 was selected for determining maximum pressure in the containment. On break of pipes in the primary and secondary loops, reactor coolant and steam in the secondary loop would discharge into the containment, rapidly and simultaneously. The pressure and temperature would reach the utmost level in the containment. Such an occurrence is not considered in the design, and the accident is beyond the design basis and irreparable damage to the containment is possible.

Given the importance of analysis of containment parameters, several studies were done to simulate LOCA accidents and find the reactor containment response. In some studies, the containment analysis program GOTHIC was used to simulate containment [1-3]. L. C. Dai et al. [4] analyzed MSLB for large dry containment using the GOTHIC code and performing the dry containment pressure and temperature analysis. The blowdown mass and energy data of the MSLB, which are tabulated in the Final Safety Analysis Report (FSAR), was used as the boundary conditions. The calculated containment pressure and temperature in the accident were compared with the FSAR results. Only two thermal hydraulic nodes were considered in their work and the condensation heat transfer between the vapor and wall surface was calculated by a simple empirical correlation. According to the most common empirical correlations, the heat transfer coefficient in a blowdown is a linear function of time, from zero to the value given by the Tagami correlation, and after the blowdown, heat transfer coefficient is given by the Uchida correlation [5].

Noori-Kalkhoran et al. [6] simulated a double ended cold leg LOCA in VVER-1000 containment as the worst case of the containment pressure increase. Thermal–hydraulic parameters of the containment were analyzed by the CONTAIN 2.0 code and a single cell model without considering condensation. Noori-Kalkhoran et al. [7] also investigated short- and long-term response of VVER-1000 containment to LOCA, considering the gap between the steel and concrete shells.

Large break LOCA in an AP1000 reactor containment was studied by Rahim et al. [8]. The containment was divided into two cells. They considered convection as a heat transfer mechanism but ignored the condensation effects. Pressurization in LOCA of an AP1000 reactor containment of three control volumes was simulated later by Rahim et al. [9].

In this paper, multicellular lumped parameter model is employed. The containment is divided into 16 cells. We use general physical models and numerical methods, which can be applied to other reactor containments.

The main assumptions in the simulation are: 1) the whole containment volume is divided into the analyzed cells; 2) each analyzed cell consist of two-component, two-phase mixture of liquid water, water-vapor, and noncondensable gases; 3) Mass and energy exchange between the cells is performed via connections; 4) the ideal gas law is used for non-condensing gases. The specific heat and other characteristics of gases are given as temperature functions. Unlike non-condensing gases, steam is considered as a real gas in saturation or superheated state, i.e. its density and enthalpy are functions of temperature and pressure; and 5) there are heat structures in some analyzed cells. Heat and mass transfer (by means of convection and condensation) between the cells and these heat structures is considered. The released mass and energy due to LOCA and MSLB are directly adopted from the FSAR analyses and are considered as initial data. In addition, the atmosphere condition outside the containment steel shell is assumed as a boundary condition for simulation zone.

2 Containment model

The VVER1000/V446 reactor containment is spherical and double-layered, an inner layer of an integrated steel shell surrounded by a concrete shell. Sixteen containment cells are considered (Fig. 1). Each has two-component of two-phase mixture of liquid water, vapor and noncondensable gases. The cells may differ in temperature and may experience heat and mass transfer. Leakage into or out of the cell is allowed by inter-cell junctions. The cell specifications are given in Table 1. Approximate room positions are highlighted in yellow on the view of the containment cross-section in Fig. 1. Connection paths between cells are indicated by red arrows. Mass transfer between cells is made possible through the paths.

Fig. 1Full-screen View

(Color online) VVER1000/V446 containment nodalization

2.3 Mass and energy transfer

In the accident, mass and energy are transferred between two cells connected to each other. Water-vapor mass and energy are released from MSLB and LOCA break spot to the specific cells. Besides, condensation and convection occurring at the control volumes near the steel shell of the containment result in heat and mass removal from these cells.

2.3.1 Junction flow model

The lumped-parameter volumes are connected by junctions using a one-dimensional model for flow between volumes. Momentum is conserved by the solution of the momentum equation in the junctions connecting the volumes as follows.

For any time step, mass flow rate between all cells are calculated by[12].

$\frac{\text{d}{W}_{ij}}{\text{d}t}=\left[P_i-P_j-{\rho }_{\text{av}}g(H_j-H_i)-C_{\text{FC}}\frac{\left|W_{ij}\right|W_{ij}}{{\rho }_{\text{av}}{A}_{ij}^2}\right]\frac{A_{ij}}{L_i{}_j},$ (15)

where, C_FC and ρ_AV are irreversible flow loss coefficient and average density, respectively; L_ij and A_ij are length and cross section of flow path, respectively; and W_ij is mass flow rate between cells i and j. W_ij is limited by the chocked flow rate:

$W_{\text{cr,}ij}=(s_{ij}-s_{ji})A_{ij}{v}_{ij}{[{\gamma }_{\text{u}}{P}_{\text{u}}{\rho }_{\text{u}}{\eta }_{\text{u}}]}^{1/2}\text{ ,}$ (16)

where, ρ_u=(ρ_i+ρ_j)/2 is the flow path density, v_ij=A_contract/A_ij is generally less than unity and is defined as the ratio of the minimum area occupied by the flow streamlines to the geometric cross-sectional area of the flow path, γ_u =C_P/C_V is the ratio of specific heats in the upstream cell atmosphere; η_u=[2/(1+γ_u)]^{(γu+1)/(γu−1)} is a dimensionless parameter; and s_ij=1 when P_i−P_j >0 and s_ij =0 when P_i−P_j <0. The subscript u denotes upstream cell conditions.

2.3.2 Flashing model

The nature of the MSLB and LOCA blowdown is different. A LOCA blowdown is liquid which flashes to a mixture of liquid and saturated steam. The resulting containment condition is a mixture of saturated steam, water, and air at saturation temperature. The MSLB blowdown is superheated steam or saturated steam which expands to superheated conditions. The resulting containment condition is a mixture of superheated steam and air, and a high degree of turbulence will be present while the blowdown continues [13]. Water-vapor mass and energy released from MSLB and LOCA are obtained from reactor safety analysis report shown in Fig. 2 [10]. The injection into the containment during the first few seconds of LOCA is entirely assumed subcooled liquid water; the injection then transitions rapidly to a mixture of steam and liquid. Pressure flash model is used to calculate how much of the entering blowdown liquid is flashed into steam based on the total compartment pressure before fluid equilibrium is reached. The liquid flashing is calculated by [14]:

Fig. 2Full-screen View

(Color online) Integration mass and energy released in containment by (a) water and steam in LOCA and (b) steam in MSLB, adopted from FSAR [10]

$M_{\text{Flash}}=M_{\text{blowdown}}\left[h_{\text{blowdown}}–h_l\left(P_{\text{t}}\right)\right]/\left[h_{\text{g}}\left(P_{\text{t}}\right)-h_l\left(P_{\text{t}}\right)\right]$ (17)

where, M_Flash is mass of blowdown liquid which flashes; P_t is total pressure of cell where fracture is located, M_blowdown and h_blowdown are mass and specific enthalpy of blowdown fluid initially entering atmosphere, respectively; and h_l is specific enthalpy of saturated liquid at pressure P_t.

2.3.3. Heat transfer models

Heat transfer from the containment atmosphere occurs through a gas boundary layer and also condensation exists in presence of noncondensable gas (air) on the containment surface. Condensation occurs through a liquid layer between the surface and the boundary layer (Fig. 3). The correlations of experimental heat and mass transfer data are usually presented with dimensionless variables, which include: Re=ρ_BLLV_C/μ_BL, Gr=max(|ρ_if−ρ_b|,10⁻⁷)L³g(ρ_BL/μ_BL)²/ρ_BL, Nu=Lh_c/K_BL, Pr=μ_BLC_P,BL/K_BL and Ra=GrPr, where, μ_BL and K_BL are the viscosity and conduction heat transfer coefficient of boundary layer, respectively; D_diff,v is the diffusivity of steam; and L is the characteristic length. These dimensionless numbers will be used in the next sections.

Fig. 3Full-screen View

The influence of noncondensables on the interface resistance.

2.3.3.1 Convection heat transfer

On basis of the scaling studies, the inside of the containment shell is expected to experience a high-velocity flow of steam and air during an MSLB event and during the blowdown phase of a large LOCA event as the break jet vigorously circulates the gas. The heat and mass transfers in this period are expected to be turbulent forced or mixed convection. After the LOCA blowdown, the atmosphere is circulated less vigorously and the velocity of the steam and air flowing along the inside surface of the containment shell will be lower. This indicates that turbulent free-convection heat and mass transfer is appropriate after blowdown [15]. Natural, forced, or mixed convection heat transfer to the atmosphere is determined at a surface by the criteria given in Table 2 [16].

Table 2Full-screen View

Nusselt number correlations for convection heat transfer [16]

Regions		Convection criteria	Nusselt number
Inside the containment
	Natural	Re2<Gr, Ra>1010	Nu = 0.228 Ra0.226
	Forced	Re2>10Gr, Re>104	Nu = 0.023 Re0.8 Pr1/3
	Mixed	1.0 Gr≤Re2≤10Gr	Nu = (Re2/Gr –1)(Nuforced –Nunatural)/9 + Nunatural
Liquid film
	Wavy laminar film	Re<5800 Pr–1.06	Nu = 0.822 Re–0.22
	Turbulent films	Re≥5800 Pr–1.06	Nu = 0.0038 Ra0.4 Pr0.65
Gap & environment
	Natural	Pr≥0.7, Ra≤1011	Nu = 2+ 0.589 Ra1/4 / [1+(0.469/Pr)9/16]4/9

Show more

We assume that the liquid film is a distinct control volume with mass transfer, convection heat transfer into the free surface, and conduction to the solid surface. The Chun and Seban correlation [17] is used to model both wavy laminar and turbulent heat transfer across the film. The correlation due to Churchill [18] is selected for calculating convection and heat transfers between outer side of the concrete wall and air through the gap shown in Table 2.

2.3.3.2 Forced convection velocities for heat transfer structures

The gas velocities for finding convection Nusselt number correlations are calculated on the basis of the gas flows into and out of a cell. It is defined as the average of the inlet and outlet velocities for a cell:

$V_{\text{C}}=\left(V_{\text{in}}+V_{\text{out}}\right)/2,$ (18)

where, V_in=M_ijRT_in/(M_u A_hd p_g), with T_in being the average temperature of the inflowing gas, M_u being the molecular weight of the gas in the upstream cell, A_hd =V_g^2/3 being the hydraulic area and V_g the cell gas volume, P_g being the pressure of the downstream cell. T_in is calculated by assuming that the flow is isothermal and that all of the incoming gas streams mix with each other before they mix with the existing cell inventory. The temperature after the streams mix in cell i is given by:

$T_{\text{in}}=\frac{{\displaystyle \sum _{\text{ji}}{s}_{\text{ji}}\left|W_{\text{ji}}\right|T_{\text{u}}{C}_{\text{p,u}}}}{{\displaystyle \sum _{\text{ji}}{s}_{\text{ji}}\left|W_{\text{ji}}\right|C_{\text{p,u}}}},$ (19)

where, T_u is the upstream cell temperature and C_p,u is the upstream cell gas specific heat at constant pressure. For gas evolving from the pool surface, these are defined as the temperature and specific heat of the evolving gas, respectively.

The expression for V_out is

$V_{\text{out}}={\displaystyle \sum _{\text{ij}}\frac{s_{\text{ij}}\left|W_{\text{ij}}\right|}{A_{\text{hd}}{\rho }_{\text{g}}}}$ (20)

where ρ_g is the density of the gas in the cell. The sum in this case includes all gas outflows through flow paths [12].

2.3.3.3 Mass transfer inside containment

Condensation mass transfer is a result of a concentration gradient between a flowing steam-air gas mixture and a surface. All water that condenses goes out of the cell. The steam concentration gradient is approximated as the difference in steam partial pressure between the bulk gas and liquid surface. Condensation occurs when the bulk gas steam concentration is greater than the concentration at the liquid surface. The steam mass flux between the heat structures surface and bulk gas is obtained based Ref. [19]:

${\dot{m}}_{\text{condens}}^{″}-k_{\text{g}}{M}_{\text{V}}(P_{\text{V}}{}_{,b}-P_{\text{V,if}}),$ (21)

where, M_V is the molecular weight of vapor; P_V,b and P_V,if are the vapor pressure at the bulk and interface, respectively; k_g=Sh P_g D_{diff, V}/(RT_BL P_nm L) is the mass transfer coefficient (in mol·s⁻¹·m⁻²·Pa⁻¹) for the gas phase mass transfer, with D_diff,V=8.54×10⁻⁵T^1.82/P for T <700 K and D_diff,V=2.2×10⁻⁴T^1.675/P for T >700 K, [20] being the boundary layer vapour diffusion coefficient (in m²/s); R being the gas constant (in Pa·cm³·mol⁻¹·K⁻¹), and P_nm=(P_V,if −P_V,b)/ln[(P_g−P_V,b)/(P_g− P_V,if)] being the logarithmic mean pressure; and Sh is the Sherwood number for mass transfer, which is obtained from the correlations of Nu described in Section 2.2.3.1 by applying the heat and mass transfer analogy. The Sherwood number is the same as the Nusselt number except that for the latter the Pr number is replaced by Schmidt number Sc as Sh =Nu{Pr→Sc}.

2.3.3.4 Boundary layer properties

All properties used in the heat and mass transfer calculations must be evaluated at a boundary layer temperature:

$T_{\text{BL}}=(T_{\text{if}}+T)/2.$ (22)

When condensation from a gaseous bulk fluid is occurring at a surface, the composition of the gas boundary layer will generally differ from the bulk composition. The vapor mole fraction must be corrected for the ongoing condensation effect. The composition of the noncondensable gases in the boundary layer, however, is still assumed to be that of the bulk gas. The vapor mole fraction for the boundary layer is obtained as:

$X_{\text{V,BL}}=(X_{\text{V,if}}+X_{\text{V,b}})/2,X_{\text{V,if}}=X_{\text{V,b}}+P_{\text{sat}}(T_{\text{if}})/P_{\text{V,b}}.$ (23)

After calculation of X_V,b and X_V,if densities of interface and the boundary layer are corrected for composition. The composition-corrected densities are also corrected for the temperature using the ideal gas law. The final equations are:

${\rho }_{\text{if,X}}=\left(\frac{X_{\text{v,if}}{M}_{\text{v}}+X_{\text{nc,if}}{M}_{\text{nc}}}{M_{\text{b}}}\right)\left(\frac{T}{T_{\text{if}}}\right){\rho }_{\text{b}},{\rho }_{\text{BL,X}}=\left(\frac{X_{\text{v,BL}}{M}_{\text{v}}+X_{\text{nc,BL}}{M}_{\text{nc}}}{M_{\text{b}}}\right)\left(\frac{T}{T_{\text{BL}}}\right){\rho }_{\text{b}},$ (24)

where, M is the molecular weight.

The gas specific heat at constant pressure in the boundary layer, C_P,BL, is evaluated at the temperature T from

$C_{\text{P,BL}}=W_{\text{V,BL}}{C}_{\text{P,V}}+W_{\text{NC,BL}}{C}_{\text{P,NC}},$ (25)

where W refers to mass fractions. Mass fractions are used because heat capacities are evaluated per unit mass, not per mole. No pressure correction is applied for C_P, the pressure dependence of C_P is weaker than the temperature dependence of the other gas transport properties. When the boundary layer properties are evaluated, only the vapor viscosity and conductivity are evaluated at the temperature T_BL by calling the appropriate property routines. For the noncondensable gas mixture, a simple power-law temperature dependence is used to correct the noncondensable properties calculated for the bulk:

${\mu }_{\text{NC,BL}}={\mu }_{\text{NC,b}}{(T_{\text{BL}}/T)}^{0.75},$ (26) $k_{\text{NC,BL}}=k_{\text{NC,b}}{(T_{\text{BL}}/T)}^{0.67}.$ (27)

After corrections are made for temperature the boundary layer properties are evaluated for the boundary layer mixture of vapor and noncondensable gases by using the mixing rules [12]:

${\mu }_{\text{BL}}=\frac{X_{\text{NC,BL}}{\mu }_{\text{NC,BL}}\sqrt{M_{\text{NC}}}+X_{\text{V,BL}}{\mu }_{\text{V,BL}}\sqrt{M_{\text{V}}}}{X_{\text{NC,BL}}\sqrt{M_{\text{NC}}}+X_{\text{V,BL}}\sqrt{M_{\text{V}}}}G(X)=0,$ (28) $K_{\text{BL}}=X_{\text{NC,BL}}{K}_{\text{NC,BL}}+X_{\text{V,BL}}{K}_{\text{V,BL}}.$ (29)

3 Numerical methods

The introduced conservation equations and constitutive relationships must be solved simultaneously. The equations are extremely nonlinear, and the inserted boundary conditions (such as mass and energy released from LOCA site) happen very fast, hence the need of an appropriate numerical method to prevent divergency of solution. Fully implicit numerical solution scheme is used to solve the coupled set of conservation equations for mass, momentum, and energy in the fluid cells, together with the heat conduction equations for the thermal conductors. Conservation equations are discretized using Crank-Nicolson discretization scheme:

where, δ₁ⁿ, δ₂ⁿ, δ₃ⁿ and δ₄ⁿ indicate the right hand of Eqs. (1), (2), (3) and (15) at time t_n. Then U, p and v_WV in the above equations are substituted by corresponding expression in state equations, i.e. Eqs.(5–10). This decreases the unknown variables to masses and temperatures of the cells and mass flow rates of the junctions. They form vector X is defined as:

The complete set of governing equations for G(X)=0 is solved simultaneously for the temperature, mass and mass flow rate variables for each new time value. Some of the terms in the equation set are nonlinear in the new time variables. An approximate solution is used by applying the one-step Newton method. The method requires an initial guess X⁽⁰⁾ as input. It then computes subsequent iterates X⁽¹⁾, X⁽²⁾... that, hopefully, will converge to a solution X*of G(X)=0. The system of linear equations which is solved for each iteration is as follows:

where, the superscript k denotes the iteration number and J_G(X) is the Jacobian matrix of G(X), defined by

The resultant matrix equation can be solved for the unknown S^(k) by direct solution (Gaussian elimination/back substitution) methods and then setting X^(k+1) = X^(k)+S^(k).

4 Results and Discussion

In Fig. 4, simulation results of LOCA are compared with the VVER1000 final safety analysis report (FSAR) [10]. Due to the large number of cells (16 cells), temperature results are given in cells 2, 5, 6, 7, 10, 14, 15 and 16. The cells are selected in such a way that almost all the main reactor compartment rooms are covered. It can be seen that, in most time, the simulation results are lower than the reference result but they are converged to each other after blowdown phase of accident. This is because there are some conservative assumptions which are considered in safety analysis (for example the 2 % decrease of containment free volume) and some simplifications employed in simulation, such as the number of cells is 23 in the safety report while it is 16 in this work. As shown in Table 3, the mean difference between the simulation and FSAR temperatures is about of 3.1%, which seems acceptable regarding the conservative calculations of FSAR in the prediction of the maximum pressure value.

Fig. 4Full-screen View

(Color online) The simulated temperatures of certain cells in LOCA, compared with the results by FSAR [10]

Fig. 5 shows average pressures of the containment simulated for LOCA, compared with those of the VVER1000 FSAR. The peak pressure calculated by the presented model is about 0.37 MPa at t=20.6 s, while it is about 0.395 MPa at t=19.9 s by FSAR. There is a good agreement between the simulation and reference results, and pressures are lower than the design pressure of 0.507 MPa. As shown in Table 4, the mean difference between the two sets of pressure results is 1.5%.

Fig. 5Full-screen View

(Color online) Comparison of simulation and FSAR average pressure during LOCA

Figure 6 shows short term and long term variations of pressure and temperature, as a function of time, calculated by the presented model when LOCA and MSLB occur simultaneously. The pressure and temperature increase rapidly due to the initial blowdown event, and eventually reach their peak values of 0.645 MPa and 153.1°C at about 32 s, where they begin to decrease slowly as the internal passive heat sinks (including the containment shell) and the steam condenses. So, if LOCA and MSLB occur simultaneously, the maximum pressure exceeds the design limit (0.51 MPa). Comparing to the LOCA as the basis accident of maximum pressure calculation in safety analyses, the peak pressure and temperature increase by 0.25 MPa and 24.9 $℃$, respectively, and the time to reach the peak pressure and temperature increases from 19.9 s to 32 s. In the safety analyses of the plant, the pressure level of 0.51 MPa has been considered as the accepted safety criteria and is used in calculations of the containment peak pressure. Regarding that the maximum pressure exceeds this value, the investigated situation is a new evaluation for the containment peak pressure.

Fig. 6Full-screen View

(Color online) Short term (a) and long term (b) pressure and temperature variations during LOCA & MSLB

5 Conclusion

A general two-phase model, including condensation in the presence of noncondensable gas has been utilized to simulate containment thermodynamic behavior when LOCA and MSLB occur concurrently. An appropriate numerical method is introduced as solution to the nonlinear governing equations.

In most of the time during LOCA, simulation results are lower than the reference result but they are converged to each other after blowdown phase of accident. Predicted temperature and pressure values that are below those of the safety report can be seen in Refs.[6,7]. One reason for this discrepancy is the difference between realistic and conservative simulation. FSAR performs a conservative safety analysis and naturally calculates a higher peak pressure. The analysis of this paper is not of safety analysis type, it is an independent evaluation of highest possible thermal-hydraulic pressure in the containment, based on the best estimate approach. Other sources of error may include simplification of simulation processes, differences between thermal structures in codes, number of thermal structures and flow paths, and difference in the structure of codes.

According to containment response in terms of both LOCA and MSLB, the pressure increases rapidly due to the initial blowdown event, reaching a peak of 0.645 MPa at about 32 s, where it begins to decrease slowly due to heat removal mechanisms. The maximum pressure exceeds the design limit (0.51 MPa). In this case, the accident would be beyond the design basis and the integrity of the containment would be threatened. Whereas reactor containment is the last barrier against the radioactive material entering the environment, which makes it particularly important in maintaining its integrity, it is suggested scenario of the break in the piping of primary and secondary loops simultaneously, is selected for the maximum pressure analysis inside the containment.