1 Introduction

Reactor containment is the last barrier to prevent the radioactive materials from entering the environment. The AP1000 reactor containment behavior during the loss of coolant accident (LOCA) has been analyzed in reactor safety reports. If the PCCS acts properly, the containment will not be exposed to any risk. However, if this system does not act properly, the situation will be much worse. In this case, the accident is beyond the design basis. From the standpoint of pressure inside the containment, the worst kind of LOCA is large LOCA. A specific type of large LOCA is double ended cold leg (DECL) break, which is a total guillotine type of break in the cold leg pipe ^[1].

DECL LOCA is characterized by two distinct phases: an initial blowdown period which continues until the reactor coolant system (RCS) pressure is nearly the same as the containment pressure, and a long term cooling period when the sensible heat of the RCS and steam generators is released along with the reactor shutdown heat. This period can continue for several hours or days until the plant is recovered. The peak pressure for this event is limited by the containment volume and by energy absorption of the passive heat sink structures inside the containment, including the steel shell. The passive cooling system (PCS) water starts to become effective at removing heat from outside of the shell just about the same time as the blowdown ends ^[2].

Given the importance of containment parameter analysis, several groups simulated LOCA and MSLB accidents and the reactor containment response, using the containment analysis program of GOTHIC ^[3-5]. L.C. Dai et al. presented the results of thermal-hydraulic calculation of “Steam Line Break” analysis at the large dry containment ^[6]. In their simulation with GOTHIC code, the blowdown mass and energy data of the main steam line break tabulated in the Final Safety Analysis Report (FSAR) were used as the boundary conditions. The calculated behavior of the containment pressure and temperature of the selected case was compared with the FSAR results, and applicability of the GOTHIC dry containment model was assessed. Only two thermal hydraulic nodes were considered, and the condensation heat transfer between the vapor and wall surface was calculated by Uchida equation. Using GOTHIC 7.2, Z. Y. Hung et al. studied the responses of AP1000 containment in two accidents: a double-ended guillotine of a hot leg and a double-ended break of the main steam line ^[7]. The results were compared with the corresponding results provided in the AP1000 ‘‘Design Control Document’’ (DCD).

A study of large break LOCA in the AP1000 reactor containment was performed by F. Rahim et al. ^[1]. They divided the containment into two cells and only thermal-hydraulic analyses were taken into account. They also considered convection as a heat transfer mechanism and ignored the condensation effect. The effect of noncondensable gasses on the steam condensation over a vertical surface was investigated by J. Su et al. ^[8] The simulation was done using CFD software with the integration of user-defined subprograms’ codes to model the heat and mass transfer during condensation. The results showed good agreement between the experimental data and the model predictions.

In all these studies, only thermal-hydraulic analyses were taken into account, while the degradation of containment directly depends on its shell strength against the abrupt rise of the pressure. Therefore, for a better understanding of containment response during the LOCA accompanied with thermal-hydraulic analysis, the stress analysis of the containment shell should be taken into account, too.

In the present work, a general two-phase model for thermal-hydraulic analysis with considering all effective phenomena such as interfacial heat and mass transfer (in presence of noncondensable gases) between vapor and liquid phases were applied. Furthermore, the obtained pressure and temperature from thermal-hydraulic models were considered as boundary conditions of the ABAQUS code to achieve the temperature and stress distribution across the steel shell of the containment during the accident. A numerical method was introduced and the thermal-hydraulic governing equations were solved by the fully implicit numerical approach, an excellent method for fast transient problems such as LB-LOCA. The MATLAB code was used for the implementation of the thermal-hydraulic models.

2 Thermal hydraulic modeling

The AP1000 reactor containment is cylindrical, double-layered (Fig. 1), equipped with a passive cooling system (PCCS). The inner layer is an integrated steel shell with an elliptical ferrule and a cap surrounded by a concrete shell. Specifications of the AP1000 reactor containment and PCCS are represented in Table 1 ^[9].

Fig. 1Full-screen View

AP1000 containment.

2.1. Containment model

It is assumed that the containment has two regions, a liquid pool region containing the reactor coolant and an atmosphere region above the pool including two-component, two-phase mixture of liquid water, water-vapor, and non-condensable gasses. Each region is assumed to have a uniform temperature, but the two regions may differ in temperature.

2.1.1 Thermodynamic model

The conservation of mass and energy equation and the equations governing the thermodynamic state of the fluid inside each region are given below, and the containment temperature and pressure as a function of time can be thus obtained.

2.1.1.1 Atmosphere region

In a time step, the mass and energy of the two-phase mixture of water-vapor in the containment atmosphere region, M_A and U_A, may be modified by evaporation or condensation interactions with the liquid pool; by heat transfer through structures, or by water and vapor leakage from the rupture site into the containment. These changes in the atmosphere region can be calculated by

$d{M}_A/dt=d{M}_{\text{WV}}/dt=W_{\text{BA}}-W_{\text{P}}-W_{\text{HS}},$ (1) $\text{d}{U}_A/\text{d}t=Q_{\text{BA}}-Q_{\text{P}}-Q_{\text{HS}},$ (2)

where, M_WV is the water vapor mass of the atmosphere region, W_BA is the water vapor mass released from rupture site into the atmosphere region as determined in Section 2.1.3, W_P is the mass transfer from atmosphere to the pool, W_HS is the mass transfer to the heat structure, Q_BA is the water vapor energy released from the rupture site into the atmosphere region; Q_P is the total heat transferred from atmosphere to the pool; and Q_HS is the total heat transferred from the containment to the heat structures. The atmosphere conditions are determined as follows:

$V_{\text{A}}=M_{\text{WV}}\times v_{\text{WV}},$ (3) $M_{\text{A}}=M_{\text{WV}}+M_{\text{NC}},$ (4) $U_{\text{A}}=U_{\text{WV}}+U_{\text{NC}},$ (5)

where subscript WV is short for water vapor, v_WV is specific volume of water vapor, subscript NC denotes non-condensable gasses, V_A, M_A, U_A are volume, mass, internal energy of the atmosphere, respectively, and

$U_{\text{WV}}=M_{\text{WV}}{u}_{\text{WV}}\left(T_{\text{A}}\right),$ (6) $U_{\text{NC}}=M_{\text{NC}}{C}_{\text{v,NC}}{T}_{\text{A}},$ (7)

where T_A is the atmosphere temperature, and C_v,NC is heat capacity of non-condensable gasses at constant volume. The internal energy u_WV and enthalpy h_WV of the two-phase mixture of water-vapor in the containment atmosphere region are calculated by:

$u_{\text{WV}}=\left(1-x\right)u_{\text{f,sat}}\left(T_{\text{A}}\right)+x{u}_{\text{g,sat}}\left(T_{\text{A}}\right)，$ (8) $h_{\text{WV}}=\left(1-x\right)h_{\text{f,sat}}\left(T_{\text{A}}\right)+x{h}_{\text{g,sat}}\left(T_{\text{A}}\right)，$ (9)

where, x is the quality of two-phase mixture. Subscripts f, g and sat denote fluid, gas and saturation, respectively. u_f,sat, u_g,sat, h_f,sat and h_g,sat are obtained from the saturation curve. In this case, the total pressure is summation of partial pressures of the vapor and noncondensable (as an ideal gas).

$P=P_{\text{wv,sat}}(T_{\text{A}})+\frac{{\text{M}}_{\text{NC}}R{T}_{\text{A}}}{x{M}_{\text{WV}}{\nu }_{\text{g,sat}}(T_{\text{A}})},$ (10)

where, P is total pressure and R is universal gas constantAlso, specific volume is determined from:

${\nu }_{\text{WV}}=(1-x){\nu }_{\text{f,sat}}(T_{\text{A}})+x{\nu }_{\text{g,sat}}(T_{\text{A}}).$ (11)

These equations are based on the assumptions of the Gibbs-Dalton law for vapors, in which the air is a perfect gas and all components are at the same temperature. Eqs. (3)–(11) are solved iteratively. M_NC is a constant parameter that is equal to the initial mass of the air in the containment atmosphere. The quantities U and M_WV are given by mass and energy conservation equations. Then, T_A, P and x can be calculated using determined parameters. Once the temperature was determined, the total pressure could be calculated from Eq. (10) ^[10].

2.1.1.2 Liquid region

During a time step, the mass and energy in the liquid region, M_P and U_P, may be modified by the interactions with the vapor region or by the water leakage from the rupture site into the containment during an accident. So, we have:

$d{M}_{\text{P}}/dt=W_{\text{BP}}+W_{\text{P}},$ (12) $d{U}_{\text{P}}/dt=Q_{\text{BP}}+Q_{\text{P}}.$ (13)

Using M_P and U_P in the new time step, the temperature and specific thermodynamic properties can be obtained. The liquid region is assumed to be a single-phase liquid at the saturation conditions. The liquid conditions are determined by:

$V_{\text{P}}=M_{\text{P}}{v}_{\text{f,sat}}\left(T_{\text{P}}\right),$ (14) $U_{\text{P}}=M_{\text{P}}{u}_{\text{f,sat}}\left(T_{\text{P}}\right).$ (15)

2.1.2 Heat and mass transfer model

During the accident, the mass and energy are transferred between the pool and atmosphere region by convection and condensation or evaporation. The condensation and convection occurring near the steel shell of the containment lead to heat and mass removal from the atmosphere. Heat and mass transfer between the atmosphere or pool and heat structure surfaces are normally evaluated with the aid of correlations involving various dimensionless numbers such as the Reynolds (Re), Prandtl (Pr), Grashof (Gr), Sherwood (Sh), and Schmidt (Sc) numbers which are defined as below.

$Re={\rho }_{\text{BL}}L{V}_{\text{C}}/{\mu }_{\text{BL}},$ (16) $Pr={\mu }_{\text{BL}}{C}_{\text{p,BL}}/K_{\text{BL}},$ (17) $Gr=L^{3}g\text{max}(|{\rho }_{\text{if}}-{\rho }_{\text{b}}|,{10}^{-7}){({\rho }_{\text{BL}}/{\mu }_{\text{BL}})}^2,G{r}_{\text{conv}}=L^{3}g\beta (T_{\text{s}}-T_{\text{b}})/v^2,$ (18) $Ra=GrPr,$ (19) $Sc={\mu }_{\text{BL}}{C}_{\text{p, BL}}/D_{\text{diff,v}},$ (20)

where, L is the characteristic length for the surface, V_C is convective gas velocity across the surface, g is acceleration of gravity, β coefficient of thermal expansion, v is kinematic viscosity. μ_BL, ρ_BL and K_BL are viscosity, density and conduction heat transfer coefficient of boundary layer. D_diff,v and C_p,BL are also diffusivity of steam and the gas specific heat at constant pressure in the boundary layer, respectively.. These dimensionless numbers, in turn, depend on the physical properties of the fluid which are calculated at the boundary layer temperature which is obtained from:

$T_{\text{BL}}=\left(T_{\text{if}}+T_{\text{b}}\right)/2,$ (21)

where, T_if is equal to T_P for the pool-atmosphere surface and to T_W for the steel wall surface when the condensation from a gaseous bulk fluid is occurring at the surface. The composition of gas boundary layer is generally different from the bulk composition. The vapor mole fraction must be corrected for the ongoing condensation. The composition of noncondensable gasses in the boundary layer, however, shall be the same as the bulk gas. The vapor mole fraction for the boundary layer is:

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

All boundary layer properties must be corrected according to the boundary layer composition and temperature. After that the corrections, the boundary layer properties are evaluated for the boundary layer mixture of the vapor and noncondensable gasses with the mixing rules. These are performed in a manner which is consistent with the recommendations of Ref.[11].

The total heat fluxes transferred from the atmosphere to the pool and steel shell consist of two components, Q_ch for convective heat transport across the boundary layer and Q_cm for the heat transported by the mass flux:

$Q″=Q_{{}_{\text{ch}}}^{″}+Q_{\text{cm}}^{″}.$ (24)

The convective heat flux is defined as ^[12]

$Q_{\text{ch}}^{″}=h_{\text{c}}\left(T_{\text{A}}-T_{\text{if}}\right),$ (25)

and the convective heat transfer coefficient is

$h_{\text{c}}=Nu{k}_{\text{BL}}/L,$ (26)

where, Nu is the Nusselt number given being 0.228 Ra^0.226, 0.15(GrPr)_b^1/3, and 0.27(GrPr)_b^1/4 for interfaces between the atmosphere and steel wall, atmosphere and pool (heated surface), and atmosphere and pool (cooled surface) ^[13].

The second component of the total heat flux is accounted for the heat transported by the mass flux for condensation or evaporation. The evaporation and condensation analytical equations are identical; only the direction of mass movement changes, so that only the condensation is described. The heat transported by the mass flux for condensation is given by:

$Q_{\text{cm}}^{″}=m_{\text{cond}}{h}_{\text{V,B}},$ (27)

where, h_V,B is enthalpy of vapor in the atmosphere region and m_cond is the condensation mass flux. Condensation mass transfer is a consequence of a concentration gradient between a flowing steam-air gas mixture and a surface. The steam concentration gradient is approximated as the difference in the 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 surface of the liquid. The steam mass flux between the surface and the bulk gas can be given as ^[14]:

$m_{{}_{\text{cond}}}^{″}=k_{\text{g}}{M}_{\text{V}}\left(P_{\text{V,B}}-P_{\text{V,if}}\right).$ (28)

The mass transfer coefficient k_g for the gas phase mass transfer is defined as:

$k_{\text{g}}=\frac{Sh{P}_{\text{g}}{D}_{\text{diff,v}}}{R{T}_{\text{BL}}{P}_{\text{nm}}L},\left(\frac{\text{kg mol}}{{\text{s m}}^{\text{2}}\text{pa}}\right),$ (29)

where, D_diff,v is the boundary layer vapor diffusion coefficient (in m²/s); R is the gas constant (in Pa·m³·kg⁻¹·mol⁻¹); and P_nm=(P_v,if −P_v,b)/ln[(P_g−P_v,b)/(P_g−P_v,if)] is the logarithmic mean pressure, with the parameters being described in Fig. 2.

Fig. 2Full-screen View

Influence of noncondensables on the interface resistance.

The Sherwood number for mass transfer, Sh, is obtained from the correlations for Nu by applying the heat and mass transfer analogy. The Sh number is equal to the Nu number with a difference and it is that in the Nusselt correlations, Pr number should be replaced by Schmidt number Sc:

$Sh=Nu\{Pr\to Sc\}$ (30)

The diffusion coefficient correlation is from Wilke-Lee modification of the Hirschfelder (WL-HBS) model ^[15]. The diffusivity is given approximately by,

$D_{\text{diff,v}}=8.54\times {10}^{-5}{T}^{1.82}/P,T<700K;D_{\text{diff,v}}=2.2\times {10}^{-4}{T}^{1.675}/P,T>700K,$ (31)

when q_p'', q_hs'', m_p'' and m_p'' are determined, we will have q_hs=q_hs''A_hs, q_p=q_p''A_pool, m_hs=m_hs''A_hs and m_p=m_p''A_pool.

2.1.3 Blowdown model

Water-vapor mass and energy released from LOCA obtained from reactor safety analysis report are shown in Fig. 3 ^[16]. The injection into the containment during the first few seconds of LOCA is entirely assumed sub-cooled 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 fluid has flashed into steam based on the total compartment pressure before fluid equilibrium has been reached. The equation used to determine the amount of blowdown liquid which has flashed is ^[17]:

Fig. 3Full-screen View

Mass discharges and energy released by steam and water, into the containment during the accident.

$M_{\text{BA}}=\left[h_{\text{blowdown}}-h_l\left(P_{\text{t}}\right)\right]m_{\text{blowdown}}/\left[h_{\text{g}}\left(P_{\text{t}}\right)-h_l\left(P_{\text{t}}\right)\right],$ (32) $Q_{\text{BA}}=h_{\text{g}}\left(P_{\text{t}}\right)M_{\text{BA}},$ (33)

where, M_BA is the mass of blowdown liquid which flashes, P_t is total pressure of the atmosphere region, m_blowdow and h_blowdown are the mass and specific enthalpy of blowdown fluid initially entering atmosphere, respectively. Any remaining non-flashed liquid (M_BA) and its energy (Q_BA) is then transferred to the pool region.

$M_{\text{BP}}=m_{\text{blowdown}}-M_{\text{BA}},$ (34) $Q_{\text{BP}}=M_{\text{BP}}{h}_l\left(P_{\text{t}}\right).$ (35)

2.1.4 Initialization of containment conditions

The water vapor partial pressure in atmosphere region is determined by

$P_{\text{WV}}=\phi P_{\text{sat,s}},$ (36)

where φ is relative humidity. The initial mass of water vapor in the atmosphere region is

$M_{\text{WV}}=\phi V_{\text{V}}/V_{\text{g}}.$ (37)

The initial mass of noncondensable gasses in the containment is calculated from:

$M_{\text{NC}}=V_{\text{V}}\left(P_{\text{t}}-P_{\text{WV}}\right)/\left(T_{\text{V}}{R}_{\text{NC}}\right).$ (38)

The initial energy associated with the noncondensable gasses is ^[10]:

$U_{\text{NC}}=M_{\text{NC}}·C_{\text{V,NC}}·T_{\text{A}}.$ (39)

2.2. PCCS model

The PCCS utilized the containment shell to remove the reactor-released energy to the environment under normal or accident conditions. In normal conditions, the air entering through the inlet of the annulus is heated by the heat flux from the steel shell and transported to the upper part of the annulus due to the differences in density. In this cooling mode, which is called dry cooling, the air continuously flows into the annulus inlet and removes the energy out to the atmosphere. In addition to the air cooling, water is sprayed on the outer shell of the containment to enhance the cooling capability of the PCCS under the accident conditions. In this cooling mode, or wet cooling, the sprayed water from the PCCS pool forms water film at the outer shell of the containment. The water film is heated and evaporated to the annulus.

Numerous studies regarding the PCCS design have been conducted using a lumped parameter code such as CONTAIN ^[18], MELCOR ^[17], or a CFD code, such as COMMIX [19], W-GOTHIC, and PCCSAC ^[20]. The W-GOTHIC code is especially used as an analysis tool in AP1000 design. In this paper, the lumped parameter model is used for the dry cooling mode and a simplified CFD model for the wet cooling mode.

2.2.1 Dry cooling mode

Except for the entrance region, the air flow inside the PCCS channel is turbulent elsewhere. Among the correlations derived from the experimental data, the correlation for the turbulent natural convection channel flow of Hugot’s is the closest to the flow in the cylindrical region of PCCS annulus. The correlation is expressed as:

$h_{\text{c}}=0.108\frac{k}{s}{\left(\frac{g\beta (T_{\text{W}}-T_{\infty })s^3}{\nu \alpha }\right)}^{0.325}.$ (40)

In the dome section, geometrical effects such as chimney effect are superimposed to the natural convection of the vertical channel. For the mixed convection model, it is recommended ^[21] to use the Churchill’s correlation as follows:

$h_{\text{c}}=\max \left[{(h_{\text{forced}}^3-h_{\text{free}}^3)}^{\frac{1}{3}},0.75h_{\text{forced}},0.75h_{\text{free}}\right],$ (41)

where, h_forced is the forced convection heat transfer coefficient and h_free is the natural convection heat transfer coefficient. We utilize the Colburn’s correlation for forced convection and the McAdams’ correlation for natural convection:

$h_{\text{forced}}=0.023kR{e}^{0.8}P{r}^{1/3}/D_{\text{h}};\text{ Colburn },$ (42) $h_{\text{free}}=0.13k{[g\beta (T_{\text{W}}-T_{\infty })/v\alpha ]}^{1/3};\text{ McAdams }.$ (43)

2.2.2 Wet cooling mode

The wet cooling model was established based on a simplified model developed by G. Aiello et al. ^[22]. In this model, the heat and mass balance between the falling film and the ascending air flow, including the thermal conduction in the film and the heat transfer from film to humid air flow by convection and evaporation, are considered. The computational domain of the model consists of the channel height (H) and width (s) delimited by two vertical walls, as shown schematically in Fig. 4. The liquid film experiences a progressive reduction of its flow rate and thickness. Due to evaporative mass transfer to the humid air, as liquid film descends along the hot wall, it may be completely dried out at a height of x*. For the mass transfer coefficient, the heat/mass transfer analogy was used ^[23]. Besides, the water film thickness (δ) was calculated as Ref.[24].

Fig. 4Full-screen View

A schematic diagram of the PCCS model.

3 Stress modeling

Most low-carbon steels such as SA738 are sensitive to temperature variations. At higher temperatures, reduction in strength decreases the allowable stresses in the structure. At lower temperatures, brittleness reduces performance and lifetime of the structure. The use of such steels in their safe range increases their reliability.

Mechanical properties of the steel are related to the pressure and temperature of containment and can be examined in the accident condition. The engineering stress-strain curves of SA-738 steel at different temperatures are shown in Fig. 5.

Fig. 5Full-screen View

Stress-strain curves of SA-738 Gr. B steel

As shown in Fig. 5, an increase in the temperature leads to a reduction of allowable stress of the material. The elastic modulus is directly related to the stiffness and strength of the material. The higher the elastic modulus, the higher is the yield stress in the stress-strain curve of the material. The Hooke's law, σ =εE, was used to calculate the yield stress, where, σ, ε and E are the stress, strain and modulus of elasticity, respectively. To evaluate the thermal and compressive stress variations across the thickness and along the height of the containment wall, the containment building of AP1000 NPP was modeled using the ABAQUS code. ABAQUS is a suite of powerful engineering simulation programs based on the finite element method. It contains an extensive library of elements that can model virtually any geometry and an extensive list of material models that can simulate the behavior of most typical engineering materials.

Two-dimensional axisymmetric solid model for thermal and compressive stresses field was used to estimate the transient temperature distribution and the resulting stresses in the containment structure. Only thermal load was considered for evaluating thermal stresses. In other words, the history of nodal temperature was defined to inner and outer nodes of steel structure as a boundary condition. The history of surface pressure was also defined to inner surface of steel structure as a loading condition and therefore the results of developed code was transferred to the ABAQUS stress model.

4 Numerical methods

The introduced conservation equations and constitutive relationships must be solved simultaneously. The presented equations are extremely nonlinear and also the inserted boundary conditions (such as mass and energy released from LOCA site) happen very fast in terms of time. Hence, an appropriate numerical method is required to prevent the divergence of the solution. An optimized numerical algorithm is presented in Fig. 6.

Fig. 6Full-screen View

The optimized numerical algorithm.

In the thermal-hydraulic simulation, the initial internal energies and masses associated with the noncondensable gasses and water-vapor mixture in the atmosphere are calculated using Eqs.(37)–(39). Next, the mass and energy conservation equations Eqs. (1), (2), (11) and (12) are implicitly solved for the next timestep to obtain the temperature and pressure values of the atmosphere and the pool. For this purpose, initial presumptions for the temperature of atmosphere and the pool close to the expected temperatures are made. These first estimations are close to the last known temperature. Then an iterative method is used to update atmosphere temperature. Afterwards, the pool temperature is updated based on the new temperature of the atmosphere using Newton's method. Then, the atmosphere temperature is updated again based on new temperature of pool and this process continues until |T_Aⁱ − T_Aⁱ⁻¹|<0.005°C and |T_Pⁱ – T_Pⁱ⁻¹|<0.005°C, with the superscript indicating the iteration number.

The iterative method is shown in Fig. 6. At first, the amount of mass and energy transferred from the atmosphere to the walls and pool are calculated based on the trial temperature and corresponding variations of the atmosphere energy and mass are determined during the time step. Then, the qualities of x_EB and x_VB are specified according to Eqs.(9) and (11), respectively, and are compared. If they are not equal to each other, new trial temperature is obtained by placing x_EB instead of x_VB in Eq.(11) and calculations are repeated until x_EB equals to x_VB.

When the atmosphere temperature is determined, the pool temperature (T_P) will be adjusted through a process similar to the Newton's method.

U_P and M_P are determined from Eqs.(12) and (11), respectively. The first trial temperature is the initial temperature or the result of the last time advancement. After the first evaluation of F(T_P), the trial temperature is changed by 0.5°K in the direction indicated by the sign of F(T_P) and F(T_P) will be reevaluated. Subsequent trial temperatures are determined by^[10]:

Iterations are terminated when |F(T_P)–F(T_P,old)|/(U_P/M_P)≤0.0005. Once the temperatures are determined, the total pressure of the atmosphere can be calculated from Eq. (10).

5 Method of FE simulation

Two-dimensional axisymmetric solid model for FE simulation is used to estimate the transient temperature and stress in the structure. Because of symmetrical condition in the geometry, material, loading and boundary condition, the two-dimensional simulation is used. The steel structure is meshed using CAX4R elements, a 4-node continuum bilinear axisymmetric quadrilateral elements with reduced integration.

In order to obtain the optimum number of elements in simulation, an analysis of the sensitivity of number of elements is performed. The convergence study of the number of elements is performed. The variation of the maximum displacement versus the number of elements along the thickness and length of the structure is extracted and the optimum of number of elements is defined. In Fig. 7, the mesh pattern of steel structure with optimum number of element, 10964, is shown.

Fig. 7Full-screen View

The mesh pattern of steel structure.

For definition the mechanical boundary conditions of the steel structure in polar coordinate system, the end of cylinder is fixed in all degree of freedom and this means that in radial, hoop and rotational about axial direction are restricted. The top of ellipsoidal section are defined as axial symmetric condition. This means that in radial, axial and rotational about hoop direction are restricted. Those boundary conditions are shown in Fig. 8.

Fig. 8Full-screen View

The boundary and loading conditions.

The history of nodal temperature that extracted from present self-code is defined to inner and outer nodes of steel structure as temperature boundary conditions. The history of surface pressure that extracted from present self-code, is also defined to inner surface of steel structure and outer surface is subjected to the atmosphere pressure (Fig. 8).

In Fig. 9(a), the contour of von misses stress and position of maximum stress is shown. The maximum stress is 23.2 MPa in end of cylindrical section. But it’s important that the hoop and radial stress must be extracted and compared with allowable limits in design. Therefore, the radial and axial stress in the end section of cylinder in the end of analysis are shown in Figs. 9(b) and 9(c). Also the hoop stress in whole of structure is shown in Fig. 9(d). The 3D-contour of hoop stress in the whole of structure is shown in Fig. 9(e).

Fig. 9Full-screen View

Contours of the von Misses stress and position of maximum stress (a), and radial (b) and axial (c) stress at the end cylindrical section.

6 Results and Discussion

In Fig.10, simulation results of blowdown and long-term cooling phase are compared with the AP1000 safety reports. The obtained results are in a good agreement with those of the AP1000 safety reports. Hence, the applied numerical method and physical models have satisfactory accuracy. In the long-term phase, reference results are a little lower than the simulation results, because that there are some simplifications employed in the simulation compared with the safety report.

Fig. 10Full-screen View

Comparison of the simulated pressure with the AP1000 safety reports [16] for blowdown phase (a) and long term cooling phase (b).

As can be seen in Fig.10, the peak pressure calculated by the presented model is about 0.483 MPa and this value is about 0.481 MPa in the reference results ^[16]. It seems that most of the time there is a good agreement between the simulation and reference results and the pressures are lower than the design pressure, which is 0.507 MPa.

According to the containment response, as shown in Fig.10, the pressure increases very rapidly due to the initial blowdown event and then the pressure decreases as the internal passive heat sinks (including the containment shell) begin to condense the steam. The pressure starts to increase again as the RCS and steam generator sensible heat is released into the containment. The pressure eventually reaches a peak at about 1800 seconds after the event’s initiation. At this time, evaporation of the PCS water cools the containment shell effectively and the pressure begins to decrease slowly.

Fig. 11 shows the variations of the temperature and pressure calculated by the presented model as a function of time when the PCCS is inactivated. As can be seen, in this case, the pressure can be reached to 0.55 MPa which is more than the design pressure.

Fig. 11Full-screen View

Variations of temperature (a) and pressure (b) inside the containment during LOCA when the PCCS system is inactivated.

According to the temperature distribution along the wall thickness at different seconds (Fig. 12a), the skin temperature reaches to its highest value at 1000 s. The elastic modulus varying with temperature (Fig. 12b) is obtained using the stress-strain curve. An increase in the temperature leads to the reduction of this parameter. Elastic modulus is directly related to the strength of steel.

Fig. 12Full-screen View

Temperature distribution in the containment wall at different seconds (a) and the E-T curve for SA-738 Gr.B steel (b).

The variations in elastic modulus in accident can be investigated using the obtained parameters. The results are presented in Fig. 13. As we expected, the temperature rises from the inner to the outer surface and decreases with increasing time, which in turn decreases the Young's modulus. The Hooke's law can be used to calculate the strength (yield stress) at the minimum elastic modulus: σ= 2.04×0.0019=0.3876 GPa.

Fig. 13Full-screen View

Variations of elastic modulus during the accident across the steel shell of containment.

The hoop and radial stress variations along depth of the containment wall (from the inner to the outer wall surface) at different seconds are shown in Fig. 14. The hoop tensile stress and compressive stress are almost the same. The maximum compressive stress is about of 0.60 MPa in the inner surface and the maximum tensile stress is about of 0.60 MPa in the outer surface. The maximum radial stress is 42 KPa in the inner and the outer surface.

Fig. 14Full-screen View

Variation of the hoop stress (a) and radial stress across the inner containment wall (b) in the accident at different seconds.

Compressive stress in the inner surface is shown in Fig. 15 as a function of time. The greatest amount of compressive stress after 1000 seconds of loading is about 60 MPa which is less than the calculated yield stress.

Fig. 15Full-screen View

Compressive stress variation across the inner containment wall.

7 Conclusion

In this study, a general two-phase model including condensation and evaporation phenomena in the presence of noncondensable gas was utilized. An optimized numerical method was introduced as a solution of the nonlinear governing equations. The result of the predicted pressure was in a good agreement for blowdown phase. In long-term phase, reference results are a little lower than the simulation results. This is because of the fact that there are some simplifications employed in the simulation compared with safety report, but the containment peak pressure was predicted accurately.

According to containment response at the moment of LOCA, the pressure increases very rapidly due to the initial blowdown event. If PCCS is activated, the containment pressure eventually reaches a peak at about 1800 seconds after the initiation of the event. At this time, evaporation of the PCS water effectively cools the containment shell and the pressure begins to decrease slowly.

If the PCCS is failed, the secondary peak pressure exceeds the design value (5.07 bar) and reaches about 0.55 MPa. In this case, the accident would be beyond the design basis. However, the FE simulation showed that the stress would still be lower than the design limit stresses. The contour of von misses stress show that the maximum stress is happened in end of cylindrical section. Hence, in the critical area (end of cylinder) hoop and radial stress were extracted and compared with allowable limits in design. Therefore, in such accident, the extracted stresses would still be lower than yield stress value and no risk would threaten the integrity of the containment.