Thermal-hydraulic and stress analysis of AP1000 reactor containment during LOCA in dry cooling mode

Thermal-hydraulic and stress analysis of AP1000 reactor containment during LOCA in dry cooling mode 1 first-author 1 corresp Corresponding author. E-mail address: sa.talebi@aut.ac.ir Corresponding author. E-mail address: sa.talebi@aut.ac.ir sa.talebi@aut.ac.ir 2 1 Department of Energy Engineering and Physics, Amirkabir University of Technology (Tehran Polytechnic), 424 Hafez Avenue, P.O. Box 15875-4413, Tehran, Iran. Department of Energy Engineering and Physics, Amirkabir University of Technology (Tehran Polytechnic), 424 Hafez Avenue, P.O. Box 15875-4413, Tehran, Iran. Young Researchers and Elite Club, Islamic Azad University, South Tehran Branch, Tehran, Iran. Young Researchers and Elite Club, Islamic Azad University, South Tehran Branch, Tehran, Iran. Some kinds of break in the reactor coolant system may cause the coolant to exit rapidly from the failure site, which leads to the loss of coolant accident (LOCA). In this paper, a stress analysis of an AP1000 reactor containment is performed in an LOCA, with the passive containment cooling system (PCCS) being available and not available for cooling the wall’s containment. The variations in the mechanical properties of the wall’s containment, including elastic modulus, strength, and stress are analyzed using the ABAQUS code. A general two-phase model is applied for modeling thermal-hydraulic behavior inside the containment. Obtained pressure and temperature from thermal-hydraulic models are considered as boundary conditions of the ABAQUS code to obtain distributions of temperature and stress across steel shell of the containment in the accident. The results indicate that if the PCCS fails the peak pressure inside the containment exceeds the design value. However, the stress would still be lower than the yield stress value, and no risk would threaten the integrity of the containment. Two-phase flow LOCA Containment integrity AP1000 reactor. 1 Introduction 1 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 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 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 3 - 5 5 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 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 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 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 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 1 The AP1000 reactor containment is cylindrical, double-layered ( Fig. 1 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 Table 1 [ 9 9 ] . 1 1 Compartments 1 1 Parameters 1 1 Values 1 1 Containment 1 1 Total volume /m 3 1 1 58969.067 1 1   1 1 Thickness /m 1 1 0.0444 1 1   1 1 Internal radius /m 1 1 19.812 1 1   1 1 Material of vessel 1 1 SA738, G.B. 1 1   1 1 Active height of heat transfer /m 1 1 47.8209 1 1   1 1 Initial pressure /Pa 1 1 108.25 1 1   1 1 Initial temperature /°C 1 1 48.89 1 1   1 1 Relative humidity 1 1 0.1 1 1 PCCS annulus 1 1 Total volume / m 3 1 1 5967 1 1   1 1 Channel height /m 1 1 27.04 1 1   1 1 Channel width /m 1 1 0.2762 1 1   1 1 Channel area /m 2 1 1 40 1 1   1 1 Initial temperature /°C 1 1 46.11 1 1   1 1 Initial pressure /Pa 1 1 101.36 1 1   1 1 Wetting fraction 1 1 0.7 1 1   1 1 Inlet air mass flow rate /kg·s −1 1 1 788.05 2.1. Containment model 2 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 3 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 4 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 <math display="block" id="M1"><mrow><mi>d</mi><msub><mi>M</mi><mi>A</mi></msub><mo>/</mo><mi>d</mi><mi>t</mi><mo>=</mo><mi>d</mi><msub><mi>M</mi><mrow><mtext>WV</mtext></mrow></msub><mo>/</mo><mi>d</mi><mi>t</mi><mo>=</mo><msub><mi>W</mi><mrow><mtext>BA</mtext></mrow></msub><mo>−</mo><msub><mi>W</mi><mtext>P</mtext></msub><mo>−</mo><msub><mi>W</mi><mrow><mtext>HS</mtext></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M2"><mrow><mtext>d</mtext><msub><mi>U</mi><mi>A</mi></msub><mo>/</mo><mtext>d</mtext><mi>t</mi><mo>=</mo><msub><mi>Q</mi><mrow><mtext>BA</mtext></mrow></msub><mo>−</mo><msub><mi>Q</mi><mtext>P</mtext></msub><mo>−</mo><msub><mi>Q</mi><mrow><mtext>HS</mtext></mrow></msub><mo>,</mo></mrow></math> 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: <math display="block" id="M3"><mrow><msub><mi>V</mi><mtext>A</mtext></msub><mo>=</mo><msub><mi>M</mi><mrow><mtext>WV</mtext></mrow></msub><mo>×</mo><msub><mi>v</mi><mrow><mtext>WV</mtext></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M4"><mrow><msub><mi>M</mi><mtext>A</mtext></msub><mo>=</mo><msub><mi>M</mi><mrow><mtext>WV</mtext></mrow></msub><mo>+</mo><msub><mi>M</mi><mrow><mtext>NC</mtext></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M5"><mrow><msub><mi>U</mi><mtext>A</mtext></msub><mo>=</mo><msub><mi>U</mi><mrow><mtext>WV</mtext></mrow></msub><mo>+</mo><msub><mi>U</mi><mrow><mtext>NC</mtext></mrow></msub><mo>,</mo></mrow></math> 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 <math display="block" id="M6"><mrow><msub><mi>U</mi><mrow><mtext>WV</mtext></mrow></msub><mo>=</mo><msub><mi>M</mi><mrow><mtext>WV</mtext></mrow></msub><msub><mi>u</mi><mrow><mtext>WV</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>A</mtext></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>,</mo></mrow></math> <math display="block" id="M7"><mrow><msub><mi>U</mi><mrow><mtext>NC</mtext></mrow></msub><mo>=</mo><msub><mi>M</mi><mrow><mtext>NC</mtext></mrow></msub><msub><mi>C</mi><mrow><mtext>v,NC</mtext></mrow></msub><msub><mi>T</mi><mtext>A</mtext></msub><mo>,</mo></mrow></math> 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: <math display="block" id="M8"><mrow><msub><mi>u</mi><mrow><mtext>WV</mtext></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>x</mi></mrow><mo>)</mo></mrow><msub><mi>u</mi><mrow><mtext>f,sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>A</mtext></msub></mrow><mo>)</mo></mrow><mo>+</mo><mi>x</mi><msub><mi>u</mi><mrow><mtext>g,sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>A</mtext></msub></mrow><mo>)</mo></mrow><mi>，</mi></mrow></math> <math display="block" id="M9"><mrow><msub><mi>h</mi><mrow><mtext>WV</mtext></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>x</mi></mrow><mo>)</mo></mrow><msub><mi>h</mi><mrow><mtext>f,sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>A</mtext></msub></mrow><mo>)</mo></mrow><mo>+</mo><mi>x</mi><msub><mi>h</mi><mrow><mtext>g,sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>A</mtext></msub></mrow><mo>)</mo></mrow><mi>，</mi></mrow></math> 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). <math display="block" id="M10"><mrow><mi>P</mi><mo>=</mo><msub><mi>P</mi><mrow><mtext>wv,sat</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>A</mtext></msub><mo stretchy="false">)</mo><mo>+</mo><mfrac><mrow><msub><mtext>M</mtext><mrow><mtext>NC</mtext></mrow></msub><mi>R</mi><msub><mi>T</mi><mtext>A</mtext></msub></mrow><mrow><mi>x</mi><msub><mi>M</mi><mrow><mtext>WV</mtext></mrow></msub><msub><mi>ν</mi><mrow><mtext>g,sat</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>A</mtext></msub><mo stretchy="false">)</mo></mrow></mfrac><mo>,</mo></mrow></math> where, P is total pressure and R is universal gas constantAlso, specific volume is determined from: <math display="block" id="M11"><mrow><msub><mi>ν</mi><mrow><mtext>WV</mtext></mrow></msub><mo>=</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>ν</mi><mrow><mtext>f,sat</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>A</mtext></msub><mo stretchy="false">)</mo><mo>+</mo><mi>x</mi><msub><mi>ν</mi><mrow><mtext>g,sat</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>A</mtext></msub><mo stretchy="false">)</mo><mo>.</mo></mrow></math> 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 10 ] . 2.1.1.2 Liquid region 4 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: <math display="block" id="M12"><mrow><mi>d</mi><msub><mi>M</mi><mtext>P</mtext></msub><mo>/</mo><mi>d</mi><mi>t</mi><mo>=</mo><msub><mi>W</mi><mrow><mtext>BP</mtext></mrow></msub><mo>+</mo><msub><mi>W</mi><mtext>P</mtext></msub><mo>,</mo></mrow></math> <math display="block" id="M13"><mrow><mi>d</mi><msub><mi>U</mi><mtext>P</mtext></msub><mo>/</mo><mi>d</mi><mi>t</mi><mo>=</mo><msub><mi>Q</mi><mrow><mtext>BP</mtext></mrow></msub><mo>+</mo><msub><mi>Q</mi><mtext>P</mtext></msub><mo>.</mo></mrow></math> 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: <math display="block" id="M14"><mrow><msub><mi>V</mi><mtext>P</mtext></msub><mo>=</mo><msub><mi>M</mi><mtext>P</mtext></msub><msub><mi>v</mi><mrow><mtext>f,sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>P</mtext></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>,</mo></mrow></math> <math display="block" id="M15"><mrow><msub><mi>U</mi><mtext>P</mtext></msub><mo>=</mo><msub><mi>M</mi><mtext>P</mtext></msub><msub><mi>u</mi><mrow><mtext>f,sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>P</mtext></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>.</mo></mrow></math> 2.1.2 Heat and mass transfer model 3 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. <math display="block" id="M16"><mrow><mi>R</mi><mi>e</mi><mo>=</mo><msub><mi>ρ</mi><mrow><mtext>BL</mtext></mrow></msub><mi>L</mi><msub><mi>V</mi><mtext>C</mtext></msub><mo>/</mo><msub><mi>μ</mi><mrow><mtext>BL</mtext></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M17"><mrow><mi>P</mi><mi>r</mi><mo>=</mo><msub><mi>μ</mi><mrow><mtext>BL</mtext></mrow></msub><msub><mi>C</mi><mrow><mtext>p,BL</mtext></mrow></msub><mo>/</mo><msub><mi>K</mi><mrow><mtext>BL</mtext></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M18"><mrow><mi>G</mi><mi>r</mi><mo>=</mo><msup><mi>L</mi><mn>3</mn></msup><mi>g</mi><mtext>max</mtext><mo stretchy="false">(</mo><mo>|</mo><msub><mi>ρ</mi><mrow><mtext>if</mtext></mrow></msub><mo>−</mo><msub><mi>ρ</mi><mtext>b</mtext></msub><mo>|</mo><mo>,</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>7</mn></mrow></msup><mo stretchy="false">)</mo><msup><mrow><mo stretchy="false">(</mo><msub><mi>ρ</mi><mrow><mtext>BL</mtext></mrow></msub><mo>/</mo><msub><mi>μ</mi><mrow><mtext>BL</mtext></mrow></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo>,</mo><mi>G</mi><msub><mi>r</mi><mrow><mtext>conv</mtext></mrow></msub><mo>=</mo><msup><mi>L</mi><mn>3</mn></msup><mi>g</mi><mi>β</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>s</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>b</mtext></msub><mo stretchy="false">)</mo><mo>/</mo><msup><mi>v</mi><mn>2</mn></msup><mo>,</mo></mrow></math> <math display="block" id="M19"><mrow><mi>R</mi><mi>a</mi><mo>=</mo><mi>G</mi><mi>r</mi><mi>P</mi><mi>r</mi><mo>,</mo></mrow></math> <math display="block" id="M20"><mrow><mi>S</mi><mi>c</mi><mo>=</mo><msub><mi>μ</mi><mrow><mtext>BL</mtext></mrow></msub><msub><mi>C</mi><mrow><mtext>p, BL</mtext></mrow></msub><mo>/</mo><msub><mi>D</mi><mrow><mtext>diff,v</mtext></mrow></msub><mo>,</mo></mrow></math> 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: <math display="block" id="M21"><mrow><msub><mi>T</mi><mrow><mtext>BL</mtext></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>if</mtext></mrow></msub><mo>+</mo><msub><mi>T</mi><mtext>b</mtext></msub></mrow><mo>)</mo></mrow><mo>/</mo><mn>2</mn><mtext> </mtext><mo>,</mo></mrow></math> 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: <math display="block" id="M22"><mrow><msub><mi>X</mi><mrow><mtext>V,BL</mtext></mrow></msub><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>X</mi><mrow><mtext>V,if</mtext></mrow></msub><mo>+</mo><msub><mi>X</mi><mrow><mtext>V,b</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>/</mo><mn>2</mn><mtext> </mtext><mo>,</mo></mrow></math> <math display="block" id="M23"><mrow><msub><mi>X</mi><mrow><mtext>V,if</mtext></mrow></msub><mo>=</mo><msub><mi>X</mi><mrow><mtext>V,b</mtext></mrow></msub><msub><mi>P</mi><mrow><mtext>sat</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mrow><mtext>if</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>/</mo><msub><mi>P</mi><mrow><mtext>V,b</mtext></mrow></msub><mo>.</mo></mrow></math> 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 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: <math display="block" id="M24"><mrow><mi>Q</mi><mo>″</mo><mo>=</mo><msubsup><mi>Q</mi><mrow><msub><mrow/><mrow><mtext>ch</mtext></mrow></msub></mrow><mo>″</mo></msubsup><mo>+</mo><msubsup><mi>Q</mi><mrow><mtext>cm</mtext></mrow><mo>″</mo></msubsup><mo>.</mo></mrow></math> The convective heat flux is defined as [ 12 12 ] <math display="block" id="M25"><mrow><msubsup><mi>Q</mi><mrow><mtext>ch</mtext></mrow><mo>″</mo></msubsup><mo>=</mo><msub><mi>h</mi><mtext>c</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>A</mtext></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>if</mtext></mrow></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>,</mo></mrow></math> and the convective heat transfer coefficient is <math display="block" id="M26"><mrow><msub><mi>h</mi><mtext>c</mtext></msub><mo>=</mo><mi>N</mi><mi>u</mi><msub><mi>k</mi><mrow><mtext>BL</mtext></mrow></msub><mo>/</mo><mi>L</mi><mo>,</mo></mrow></math> 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 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: <math display="block" id="M27"><mrow><msubsup><mi>Q</mi><mrow><mtext>cm</mtext></mrow><mo>″</mo></msubsup><mo>=</mo><msub><mi>m</mi><mrow><mtext>cond</mtext></mrow></msub><msub><mi>h</mi><mrow><mtext>V,B</mtext></mrow></msub><mo>,</mo></mrow></math> 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 14 ] : <math display="block" id="M28"><mrow><msubsup><mi>m</mi><mrow><msub><mrow/><mrow><mtext>cond</mtext></mrow></msub></mrow><mo>″</mo></msubsup><mo>=</mo><msub><mi>k</mi><mtext>g</mtext></msub><msub><mi>M</mi><mtext>V</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mrow><mtext>V,B</mtext></mrow></msub><mo>−</mo><msub><mi>P</mi><mrow><mtext>V,if</mtext></mrow></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>.</mo></mrow></math> The mass transfer coefficient k g for the gas phase mass transfer is defined as: <math display="block" id="M29"><mrow><msub><mi>k</mi><mtext>g</mtext></msub><mo>=</mo><mfrac><mrow><mi>S</mi><mi>h</mi><msub><mi>P</mi><mtext>g</mtext></msub><msub><mi>D</mi><mrow><mtext>diff,v</mtext></mrow></msub></mrow><mrow><mi>R</mi><msub><mi>T</mi><mrow><mtext>BL</mtext></mrow></msub><msub><mi>P</mi><mrow><mtext>nm</mtext></mrow></msub><mi>L</mi></mrow></mfrac><mtext> </mtext><mo>,</mo><mrow><mo>(</mo><mrow><mfrac><mrow><mtext>kg mol</mtext></mrow><mrow><msup><mrow><mtext>s m</mtext></mrow><mtext>2</mtext></msup><mtext>pa</mtext></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where, D diff,v is the boundary layer vapor diffusion coefficient (in m 2 /s); R is the gas constant (in Pa·m 3 ·kg −1 ·mol −1 ); 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. 2 . 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 : <math display="block" id="M30"><mrow><mi>S</mi><mi>h</mi><mo>=</mo><mi>N</mi><mi>u</mi><mo>{</mo><mi>P</mi><mi>r</mi><mo>→</mo><mi>S</mi><mi>c</mi><mo>}</mo></mrow></math> The diffusion coefficient correlation is from Wilke-Lee modification of the Hirschfelder (WL-HBS) model [ 15 15 ] . The diffusivity is given approximately by, <math display="block" id="M31"><mrow><msub><mi>D</mi><mrow><mtext>diff,v</mtext></mrow></msub><mo>=</mo><mn>8.54</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>5</mn></mrow></msup><msup><mi>T</mi><mrow><mn>1.82</mn></mrow></msup><mo>/</mo><mi>P</mi><mo>,</mo><mi>T</mi><mo><</mo><mn>700</mn><mtext> </mtext><mi>K</mi><mo>;</mo><mtext> </mtext><msub><mi>D</mi><mrow><mtext>diff,v</mtext></mrow></msub><mo>=</mo><mn>2.2</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>4</mn></mrow></msup><msup><mi>T</mi><mrow><mn>1.675</mn></mrow></msup><mo>/</mo><mi>P</mi><mo>,</mo><mi>T</mi><mo>></mo><mn>700</mn><mtext> </mtext><mi>K</mi><mo>,</mo></mrow></math> 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 3 Water-vapor mass and energy released from LOCA obtained from reactor safety analysis report are shown in Fig. 3 Fig. 3 [ 16 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 17 ] : <math display="block" id="M32"><mrow><msub><mi>M</mi><mrow><mtext>BA</mtext></mrow></msub><mo>=</mo><mtext> </mtext><mrow><mo>[</mo><mrow><msub><mi>h</mi><mrow><mtext>blowdown</mtext></mrow></msub><mo>−</mo><msub><mi>h</mi><mi>l</mi></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><msub><mi>m</mi><mrow><mtext>blowdown</mtext></mrow></msub><mo>/</mo><mrow><mo>[</mo><mrow><msub><mi>h</mi><mtext>g</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow><mo>−</mo><msub><mi>h</mi><mi>l</mi></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mtext> </mtext><mo>,</mo></mrow></math> <math display="block" id="M33"><mrow><msub><mi>Q</mi><mrow><mtext>BA</mtext></mrow></msub><mo>=</mo><msub><mi>h</mi><mtext>g</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow><msub><mi>M</mi><mrow><mtext>BA</mtext></mrow></msub><mo>,</mo></mrow></math> 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. <math display="block" id="M34"><mrow><msub><mi>M</mi><mrow><mtext>BP</mtext></mrow></msub><mo>=</mo><msub><mi>m</mi><mrow><mtext>blowdown</mtext></mrow></msub><mo>−</mo><msub><mi>M</mi><mrow><mtext>BA</mtext></mrow></msub><mo>,</mo></mrow></math> <math display="block" id="M35"><mrow><msub><mi>Q</mi><mrow><mtext>BP</mtext></mrow></msub><mo>=</mo><msub><mi>M</mi><mrow><mtext>BP</mtext></mrow></msub><msub><mi>h</mi><mi>l</mi></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>.</mo></mrow></math> 2.1.4 Initialization of containment conditions 3 The water vapor partial pressure in atmosphere region is determined by <math display="block" id="M36"><mrow><msub><mi>P</mi><mrow><mtext>WV</mtext></mrow></msub><mo>=</mo><mi>φ</mi><msub><mi>P</mi><mrow><mtext>sat,s</mtext></mrow></msub><mo>,</mo></mrow></math> where φ is relative humidity. The initial mass of water vapor in the atmosphere region is <math display="block" id="M37"><mrow><msub><mi>M</mi><mrow><mtext>WV</mtext></mrow></msub><mo>=</mo><mi>φ</mi><msub><mi>V</mi><mtext>V</mtext></msub><mo>/</mo><msub><mi>V</mi><mtext>g</mtext></msub><mo>.</mo></mrow></math> The initial mass of noncondensable gasses in the containment is calculated from: <math display="block" id="M38"><mrow><msub><mi>M</mi><mrow><mtext>NC</mtext></mrow></msub><mo>=</mo><msub><mi>V</mi><mtext>V</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>P</mi><mtext>t</mtext></msub><mo>−</mo><msub><mi>P</mi><mrow><mtext>WV</mtext></mrow></msub></mrow><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mrow><msub><mi>T</mi><mtext>V</mtext></msub><msub><mi>R</mi><mrow><mtext>NC</mtext></mrow></msub></mrow><mo>)</mo></mrow><mtext> </mtext><mo>.</mo></mrow></math> The initial energy associated with the noncondensable gasses is [ 10 10 ] : <math display="block" id="M39"><mrow><msub><mi>U</mi><mrow><mtext>NC</mtext></mrow></msub><mo>=</mo><msub><mi>M</mi><mrow><mtext>NC</mtext></mrow></msub><mo>·</mo><msub><mi>C</mi><mrow><mtext>V,NC</mtext></mrow></msub><mo>·</mo><msub><mi>T</mi><mtext>A</mtext></msub><mo>.</mo></mrow></math> 2.2. PCCS model 2 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 18 ] , MELCOR [ 17 17 ] , or a CFD code, such as COMMIX [ 19 19 ], W-GOTHIC, and PCCSAC [ 20 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 3 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: <math display="block" id="M40"><mrow><msub><mi>h</mi><mtext>c</mtext></msub><mo>=</mo><mn>0.108</mn><mfrac><mi>k</mi><mi>s</mi></mfrac><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mi>g</mi><mi>β</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>W</mtext></msub><mo>−</mo><msub><mi>T</mi><mi>∞</mi></msub><mo stretchy="false">)</mo><msup><mi>s</mi><mn>3</mn></msup></mrow><mrow><mi>ν</mi><mi>α</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>0.325</mn></mrow></msup><mo>.</mo></mrow></math> 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 21 ] to use the Churchill’s correlation as follows: <math display="block" id="M41"><mrow><msub><mi>h</mi><mtext>c</mtext></msub><mo>=</mo><mi>max</mi><mrow><mo>[</mo><mrow><msup><mrow><mo stretchy="false">(</mo><msubsup><mi>h</mi><mrow><mtext>forced</mtext></mrow><mn>3</mn></msubsup><mo>−</mo><msubsup><mi>h</mi><mrow><mtext>free</mtext></mrow><mn>3</mn></msubsup><mo stretchy="false">)</mo></mrow><mrow><mfrac><mn>1</mn><mn>3</mn></mfrac></mrow></msup><mo>,</mo><mn>0.75</mn><mtext> </mtext><msub><mi>h</mi><mrow><mtext>forced</mtext></mrow></msub><mo>,</mo><mn>0.75</mn><mtext> </mtext><msub><mi>h</mi><mrow><mtext>free</mtext></mrow></msub></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> 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: <math display="block" id="M42"><mrow><msub><mi>h</mi><mrow><mtext>forced</mtext></mrow></msub><mo>=</mo><mn>0.023</mn><mtext> </mtext><mi>k</mi><mtext> </mtext><mi>R</mi><msup><mi>e</mi><mrow><mn>0.8</mn></mrow></msup><mi>P</mi><msup><mi>r</mi><mrow><mn>1</mn><mo>/</mo><mn>3</mn></mrow></msup><mo>/</mo><msub><mi>D</mi><mtext>h</mtext></msub><mo>;</mo><mtext> Colburn </mtext><mo>,</mo></mrow></math> <math display="block" id="M43"><mrow><msub><mi>h</mi><mrow><mtext>free</mtext></mrow></msub><mo>=</mo><mn>0.13</mn><mi>k</mi><msup><mrow><mo stretchy="false">[</mo><mi>g</mi><mi>β</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>W</mtext></msub><mo>−</mo><msub><mi>T</mi><mi>∞</mi></msub><mo stretchy="false">)</mo><mo>/</mo><mi>v</mi><mi>α</mi><mo stretchy="false">]</mo></mrow><mrow><mn>1</mn><mo>/</mo><mn>3</mn></mrow></msup><mo>;</mo><mtext> McAdams </mtext><mo>.</mo></mrow></math> 2.2.2 Wet cooling mode 3 The wet cooling model was established based on a simplified model developed by G. Aiello et al . [ 22 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 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 23 ] . Besides, the water film thickness ( δ ) was calculated as Ref.[ 24 24 ]. 3 Stress modeling 1 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. 5 . As shown in Fig. 5 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 1 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. 6 . 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 i − T A i −1 | < 0.005°C and | T P i – T P i −1 | < 0.005°C, with the superscript indicating the iteration number. The iterative method is shown in Fig. 6 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. <math display="block" id="M44"><mrow><mi>F</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>P</mtext></msub><mo stretchy="false">)</mo><mo>=</mo><msub><mi>U</mi><mtext>P</mtext></msub><mo>/</mo><msub><mi>M</mi><mtext>P</mtext></msub><mo>−</mo><msub><mi>u</mi><mrow><mtext>sat,f</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>P</mtext></msub><mo stretchy="false">)</mo><mo>.</mo></mrow></math> 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 10 ] : <math display="block" id="M45"><mrow><msub><mi>T</mi><mrow><mtext>p,new</mtext></mrow></msub><mo>=</mo><msub><mi>T</mi><mtext>p</mtext></msub><mo>−</mo><mfrac><mrow><mi>F</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>P</mtext></msub><mo stretchy="false">)</mo></mrow><mrow><mrow><mo>[</mo><mrow><mi>F</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>P</mtext></msub><mo stretchy="false">)</mo><mo>−</mo><mi>F</mi><mo stretchy="false">(</mo><msub><mi>T</mi><mrow><mtext>P,old</mtext></mrow></msub><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow><mo>/</mo><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>P</mtext></msub><mo>−</mo><msub><mi>T</mi><mrow><mtext>P,old</mtext></mrow></msub><mo stretchy="false">)</mo></mrow></mfrac><mo>.</mo></mrow></math> 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 1 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 Fig. 7 , the mesh pattern of steel structure with optimum number of element, 10964, is shown. 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. 8 . 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 Fig. 8 ). In Fig. 9(a) 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) Figs. 9(b) and 9(c) . Also the hoop stress in whole of structure is shown in Fig. 9(d) Fig. 9(d) . The 3D-contour of hoop stress in the whole of structure is shown in Fig. 9(e) Fig. 9(e) . 6 Results and Discussion 1 In Fig.10 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. As can be seen in Fig.10 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 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 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 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. According to the temperature distribution along the wall thickness at different seconds ( Fig. 12a Fig. 12a ), the skin temperature reaches to its highest value at 1000 s. The elastic modulus varying with temperature ( Fig. 12b 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. The variations in elastic modulus in accident can be investigated using the obtained parameters. The results are presented in Fig. 13 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. 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. 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. Compressive stress in the inner surface is shown in Fig. 15 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. 7 Conclusion 1 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. References F. Rahim , M. Rahgoshay , S. Mousavian , A Study of large break LOCA in the AP1000 reactor containment . Prog. Nucl. Energy . 54 , 132 - 137 ( 2012 ). DOI: 10.1016/j.pnucene.2011.07.004 http://doi.org/10.1016/j.pnucene.2011.07.004 A Study of large break LOCA in the AP1000 reactor containment NRC , Evaluation of the effect of the AP1000 enhanced shield building design on the containment response and safety analyses . Westinghouse Electric Company LLC , Pittsburgh ( 2011 ). Evaluation of the effect of the AP1000 enhanced shield building design on the containment response and safety analyses D. Papini , D. Grgić , A. Cammi , M.E. Ricotti , Analysis of different containment models for IRIS small break LOCA, using GOTHIC and RELAP5 codes . Nucl. Eng. Des . 241 , 1152 - 1164 ( 2011 ). DOI: 10.1016/j.nucengdes.2010.06.016 http://doi.org/10.1016/j.nucengdes.2010.06.016 Analysis of different containment models for IRIS small break LOCA, using GOTHIC and RELAP5 codes Y. Chen , Y. Yuann , L. Dai , Lungmen ABWR containment analyses during short-term main steam line break LOCA using GOTHIC . Nucl. Eng. Des . 247 , 106 - 115 ( 2012 ). DOI: 10.1016/j.nucengdes.2012.02.012 http://doi.org/10.1016/j.nucengdes.2012.02.012 Lungmen ABWR containment analyses during short-term main steam line break LOCA using GOTHIC T. Kim , J. Park , A containment analysis for SBLOCA in the refurbished Wolsong-1 nuclear power plant . Nucl. Eng. Des . 241 , 3804 - 3811 ( 2011 ). DOI: 10.1016/j.nucengdes.2011.06.040 http://doi.org/10.1016/j.nucengdes.2011.06.040 A containment analysis for SBLOCA in the refurbished Wolsong-1 nuclear power plant L.C. Dai , Y.S. Chen , Y. Yuann , Short-term pressure and temperature MSLB response analyses for large dry containment of the Maanshan nuclear power station . Nucl. Eng. Des . 280 , 86 - 93 ( 2014 ). DOI: 10.1016/j.nucengdes.2014.09.007 http://doi.org/10.1016/j.nucengdes.2014.09.007 Short-term pressure and temperature MSLB response analyses for large dry containment of the Maanshan nuclear power station Z.Y. Hung , Y.M. Ferng , W.S. Hsu et al ., Analysis of AP1000 containment passive cooling system during a loss-of-coolant accident . Ann. Nucl. Energy . 85 , 717 - 724 ( 2015 ). DOI: 10.1016/j.anucene.2015.06.027 http://doi.org/10.1016/j.anucene.2015.06.027 Analysis of AP1000 containment passive cooling system during a loss-of-coolant accident J. Su , Z. Sun , D. Zhang . Numerical analysis of steam condensation over a vertical surface in presence of air . Ann. Nucl. Energy . 72 , 268 - 276 ( 2014 ). DOI: 10.1016/j.anucene.2014.05.019 http://doi.org/10.1016/j.anucene.2014.05.019 Numerical analysis of steam condensation over a vertical surface in presence of air NRC , Final Safety Evaluation Report Related to Certification of the AP1000 Standard Design . U.S. Nuclear Regulatory Commission , Washington ( 2004 ). Final Safety Evaluation Report Related to Certification of the AP1000 Standard Design D.O.N.W. Hargroves , L.J. Metcalfe , L.L. Wheat et al ., CONTEMPT-LT/028- A computer program for predicting containment pressure- temperature response to a loss-of-coolant accident . Idaho National Engineering Laboratory , Idaho ( 1979 ). CONTEMPT-LT/028- A computer program for predicting containment pressure- temperature response to a loss-of-coolant accident K. Murata , D. Williams , J. Tills , Code manual for CONTAIN 2.0 a computer code for nuclear reactor containment analysis . Sandia National Labratories , Albuquerque, NM ( 1997 ). Code manual for CONTAIN 2.0 a computer code for nuclear reactor containment analysis S. Talebi , F. Abbasi , H. Davilu , A 2D numerical simulation of sub-cooled flow boiling at low-pressure and low-flow rates . Nucl. Eng. Des . 239 , 140 - 146 ( 2009 ). DOI: 10.1016/j.nucengdes.2008.09.007 http://doi.org/10.1016/j.nucengdes.2008.09.007 A 2D numerical simulation of sub-cooled flow boiling at low-pressure and low-flow rates Y.A. Cengel , A.J. Ghajar , Heat and Mass Transfer: Fundamentals and Applications . New York , McGraw-Hill ( 2010 ). Heat and Mass Transfer: Fundamentals and Applications F. Kreith , Principles of Heat Transfer . Harper & Row , New York ( 1973 ). Principles of Heat Transfer R. Perry , D. Green , J. Maloney , Perry’s Chemical Engineer's Handbook . McGraw-Hil , New York ( 1997 ). Perry’s Chemical Engineer's Handbook NRC , Westinghouse AP1000 Design Control Document. Rev. 19, Section 3.8 , Nuclear Regulatory Commission , Washington ( 2011 ). Westinghouse AP1000 Design Control Document. Rev. 19, Section 3.8 J. Tills , A. Notafrancesco , J. Phillips , Application of the MELCOR code to design basis PWR large dry containment analysis . Sandia National Laboratories , New Mexico ( 2009 ). Application of the MELCOR code to design basis PWR large dry containment analysis R. Vijaykumar , M. Khatib-Rahbar , Applicability of the CONTAIN code heat and mass-transfer models to asymmetrically heated vertical channels . Nucl. Technol . 128 , 313 - 326 ( 1999 ). Applicability of the CONTAIN code heat and mass-transfer models to asymmetrically heated vertical channels W.T. Sha , T.H. Chien , J.G. Sun , et al ., Analysis of large-scale tests for AP-600 passive containment cooling system . Nucl. Eng. Des . 232 , 197 - 216 ( 2004 ). DOI: 10.1016/j.nucengdes.2004.05.003 http://doi.org/10.1016/j.nucengdes.2004.05.003 Analysis of large-scale tests for AP-600 passive containment cooling system Y. Jiyang , J. Baoshan , PCCSAC: A three-dimensional code for AC600 passive containment cooling system analysis . Nucl. Sci. Eng . 142 , 230 - 236 ( 2002 ). PCCSAC: A three-dimensional code for AC600 passive containment cooling system analysis Y.D. Hwang , B.D. Chung , B.H. Cho et al ., PCCS Analysis model for the passively cooled steel containment . Nucl. Eng. Technol . 30 , 26 - 39 ( 1998 ). PCCS Analysis model for the passively cooled steel containment G. Aiello , M. Ciofalo , Natural convection cooling of a hot vertical wall wet by a falling liquid film . Int. J. Heat Mass Transf . 52 , 5954 - 5961 ( 2009 ). DOI: 10.1016/j.ijheatmasstransfer.2009.08.007 http://doi.org/10.1016/j.ijheatmasstransfer.2009.08.007 Natural convection cooling of a hot vertical wall wet by a falling liquid film C.H. Song , D.Y. Lee , S.T. Ro , Cooling enhancement in an air-cooled finned heat exchanger by thin water film . Int. J. Heat Mass Transf . 46 , 1241 - 1249 ( 2003 ). DOI: 10.1016/S0017-9310(02)00405-2 http://doi.org/10.1016/S0017-9310(02)00405-2 Cooling enhancement in an air-cooled finned heat exchanger by thin water film X. Huang , X. Cheng , Modification and application of water film model in COCOSYS for PWR’s passive containment cooling . Nucl. Eng. Des . 280 , 251 - 261 ( 2014 ). DOI: 10.1016/j.nucengdes.2014.08.026 http://doi.org/10.1016/j.nucengdes.2014.08.026 Modification and application of water film model in COCOSYS for PWR’s passive containment cooling 10.1007/s41365-017-0233-8 2016-06-18 2016-10-07 2016-10-09 2017-04-28 2017-06-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0233-8 http://dx.doi.org/10.1007/s41365-017-0233-8 . 2017 Nuclear Science and Techniques 06 28