CC BY-NC-ND
0Annals of Nuclear Energy xxx (2015) xxx-xxx
Contents lists available at ScienceDirect
Annals of Nuclear Energy
journal homepage: www.elsevier.com/locate/anucene
A two dimensional approach for temperature distribution in reactor lower head during severe accident
Zhen Cao, Xiaojing Liu *, Xu Cheng
School of Nuclear Science and Engineering, Shanghai Jiao Tong University, 800 Dong Chuan Road, Shanghai 200240, China
ARTICLE INFO
ABSTRACT
Article history:
Received 30 October 2014
Received in revised form 25 January 2015
Accepted 10 April 2015
Available online xxxx
Keywords:
Two-dimensional
Code development
In-vessel retention
Post-CHF
Vessel conduction
Boiling crisis spreading
In order to evaluate the safety margin during a postulated severe accident, a module named ASAP-2D (Accident Simulation on Pressure vessel-2 Dimensional), which can be implemented into the severe accident simulation codes (such as ATHLET-CD), is developed in Shanghai Jiao Tong University. Based on two-dimensional spherical coordinates, heat conduction equation for transient state is solved implicitly. Together with solid vessel thickness, heat flux distribution and heat transfer coefficient at outer vessel surface are obtained. Heat transfer regime when critical heat flux has been exceeded (POST-CHF regime) could be simulated in the code, and the transition behavior of boiling crisis (from spatial and temporal points of view) can be predicted.
The module is verified against a one-dimensional analytical solution with uniform heat flux distribution, and afterwards this module is applied to the benchmark illustrated in NUREG/CR-6849. Benchmark calculation indicates that maximum heat flux at outer surface of RPV could be around 20% lower than that of at inner surface due to two-dimensional heat conduction. Then a preliminary analysis is performed on the integrity of the reactor vessel for which the geometric parameters and boundary conditions are derived from a large scale advanced pressurized water reactor. Results indicate that heat flux remains lower than critical heat flux. Sensitivity analysis indicates that outer heat flux distribution is more sensitive to input heat flux distribution and the transition boiling correlation than mass flow rate in external reactor vessel cooling (ERVC) channel, and the correlation for molten vessel and ERVC coolant inlet temperature.
According to the results achieved, the new developed module shows good applicability to simulate the pressure vessel behavior during melt pool formation. Thus it can be applied for the future study of the severe accidents relating to lower head integrity.
© 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CCBY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4XI/).
1. Introduction
Fukushima accident has arisen international attention on the nuclear energy safety, especially under severe accident (SA) condition (Wittneben, 2012). Therefore in spite of core melt frequency (CMF) per reactor year for currently operating reactors has been reduced as low as 1 x 10~4 (Schulz, 2006), plenty of work needs to be done to prevent the release of the radioactive material to environment under severe accident. In order to ensure radioactive material is under control different severe accident management (SAM) strategies are developed. One of candidate strategies is in-vessel retention (IVR) strategy, which has been approved to be part of SAM for Loviisa plant (Theofanous et al., 1996). Then IVR strategy was proposed for AP600, AP1000 (Westinghouse design),
* Corresponding author. Tel.: +86 21 34207121. E-mail address: xiaojingliu@sjtu.edu.cn (X. Liu).
SWR-1000 (KERENA) and Korea APR-1400 design. In-vessel retention is an effective severe accident management strategy provided that the margin between critical heat flux and lower head heat flux is large enough.
The success of IVR strategy strongly relies on heat removal capability through external reactor vessel cooling (ERVC). In order to investigate integrity of lower head, several improvements have been done to current system codes to simulate thermal hydraulic behavior of molten materials in lower head, including MELCOR (Gauntt, 2005), SCDAP/RELAP5 (The SCDAP/RELAP5-3D Code Development Team, 2003), MAAP (MAAP code development Team, 1994), ATHLET-CD (Austregesilo, 2012). Alternatively, some codes are specially developed to evaluate integrity of the reactor pressure vessel (RPV) wall cooled from outside, such as VESSEL, MVITA (Sehgal et al., 2003), ERI-IVRAM (Esmaili and Khatib-Rahbar, 2004) and VESTA (Rempe et al., 1997). Since the temperature field of lower head is decisive factor in identifying the effectiveness of IVR
http://dx.doi.org/10.1016/j.anucene.2015.04.042 0306-4549/© 2015 The Authors. Published by Elsevier Ltd.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
2 Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
Nomenclature
General symbolsAb
inner lower head surface area (m2)
Aout outer lower head surface area (m2)
c specific heat capacity (vessel wall) (J/kg K)
cp specific heat capacity (fluid) (J/kg K)
Gj external vessel mass flow rate (kg/s)
g gravity acceleration (m/s2)
hc heat transfer coefficient for convection (W/m2 K)
hf fluid enthalpy in external vessel channel (J/kg)
hfb heat transfer coefficient for film boiling (W/m2 K)
his saturation enthalpy of liquid phase (J/kg)
hv latent heat of vaporization (J/kg)
hiv effective latent heat of vaporization (J/kg)
hnb heat transfer coefficient for nucleate boiling (J/kg)
honb enthalpy corresponding to onset nucleate boiling (J/kg)
hout heat transfer coefficient at outer vessel (W/m2 K)
hsat saturated liquid enthalpy (J/kg)
hvs saturation enthalpy of vapor phase (J/kg)
j control volume number
L characteristic length (m)
N iteration step
Nonb control volume number of first nucleate boiling volume
Nu Nusselt number
p pressure (Pa)
Pr Prandtl number
Qa heat flux for convection caused bubble agitation (W/m2)
qchf critical heat flux (W/m2)
Qe heat flux for latent heat for vapor formation (W/m2)
Qin inner vessel heat flux (W/m2)
qout outer vessel heat flux (W/m2)
Qsp heat flux for single-phase convection (W/m2)
r distance from inner vessel wall (m)
Ra Rayleigh number
Re Reynolds number
rin inner radius of lower head (m)
rout outer radius of lower head (m)
t time step (s)
T temperature (K)
Tchf wall temperature corresponding to qchf (K)
Tf fluid temperature in external vessel channel (K)
Tmelt melting temperature for vessel wall (K)
Tmfb minimum film boiling temperature (K)
Tonb wall temperature corresponding to onset nucleate
boiling (K)
Tsat saturated fluid temperature (K)
xqm mass quality
Greek symbols
a thermal diffusivity (m2/s)
b thermal expansion coefficient (K_1)
e pumping factor
h angle from bottom center of lower head (°)
k thermal conductivity (W/m K)
l dynamic viscosity (Pa s)
q density (kg/m3)
r surface tension (N/m)
Subscripts
CHF critical heat flux
in inner surface of vessel
l liquid phase
out outer surface of vessel
onb onset nucleate boiling
s saturate sate
v vapor phase
strategy, it is important to obtain relative accurate temperature field of lower head. Despite the progress that has been made in these years, there are still some parts need to be improved. Thus key phenomenon in lower head during severe accident should be considered. After an investigation of the IVR process, the important features in the IVR analysis are summarized as followed:
(1) Lower head is part of spherical shell, as is shown in Fig. 1, heat conduction in lower head occurs both r direction and h direction, besides inner surface area is smaller than outer one.
(2) The thermal conductivity of lower head material varies with
temperature.
Fig. 1. Two dimensional heat conduction.
(3) Heat transfer could be enhanced when vessel wall is partly molten due to the convection.
(4) Wide range of external vessel heat transfer correlations including POST-CHF correlations are needed to predict the surface temperature after occurrence of CHF.
Heat transfer models applied in different codes are summarized in Table 1. It can be seen that current codes have challenges in describing some of above features in lower head. To be specific, in ATHLET-CD (Austregesilo, 2012), ERI-IVRAM (Esmaili and Khatib-Rahbar, 2004) and VESTA (Rempe et al., 1997), one-dimensional heat conduction model in lower head is applied. This simplified model used to be reasonable due to their conservative estimation. However the margin to CHF is small when relatively larger thermal power reactor is applied. It is important to predict local heat flux more accurately by applying at least two-dimensional heat conduction model in lower head. In fact heat flux at outer vessel wall would be flattened due to heat conduction in h direction. Due to larger surface area, average outer surface heat flux (qout) is lower than that of inner surface (qin). This difference between qin and qout will increase with increasing lower head radius. For instance in large-scale advanced PWR in China, the difference is about 14%. However most of these codes neglect this difference. So when IVR analysis is performed, two-dimensional heat conduction model in spherical coordinate should be considered for some cases.
Apparently, thermal conductivity varies with temperature, however, in ERI-IVRAM (Esmaili and Khatib-Rahbar, 2004) and VESTA (Rempe et al., 1997) constant lower head thermal conductivity is used. Heat flux distribution in lower head would not show much
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
Table 1
Comparison of different codes.
2D heat conduction model Thermal conductivity for molten vessel Detailed external heat transfer mode
ATHLET-CD X X X
SCDAP/RELAP P X -
MELCOR P - -
MAAP P - -
VESSEL P - -
ERI-IVRAM X X X
VESTA X X X
x: has not been considered; -: model is relatively simple; model is suitable
difference if constant thermal conductivity in one-dimensional heat conduction model is used. However the difference is notable when taking two-dimensional heat conduction into consideration. This is because conductivity for colder part (like bottom of lower head) is larger than that in hotter part (like position adjacent to 90°). This effect could reduce ''focusing effect'' due to heat conduction in Q direction. Thus, it is necessary and important to consider the variation of thermal conductivity with temperature during analyzing IVR.
In addition, heat transfer could be enhanced when RPV wall is partially melted. In ATHLET-CD (Austregesilo, 2012), SCDAP/ RELAP (The SCDAP/RELAP5-3D Code Development Team, 2003) and MELCOR (Gauntt, 2005), the thickness of vessel wall would keep constant even when wall temperature exceeds melting temperature. Thermal conductivity for molten vessel is increased to some extent. However these codes have challenges in obtaining rational inner wall temperature due to lack of reasonable physical model. Thus an improved heat conduction model for molten vessel suitable for two-dimensional heat conduction is needed.
Since heat capacity of lower head is relatively large, it is possible that lower head can keep integrity even shortly after boiling crisis occurs, provided that transverse heat conduction (in Q direction in Fig. 1) in transition boiling region is sufficient. Thus suitable nucleate boiling heat transfer correlations as well as POST-CHF correlations should be considered. However, in ERI-IVRAM, only nucleate heat transfer correlation is utilized. Although MELCOR and SCDAP/RELAP employ nucleate boiling correlation, transition boiling correlation and film boiling correlation, the difference between implemented transition boiling correlations of these codes is not negligible. Thus, further study concerning POST-CHF heat transfer correlations need to be analyzed.
In order to obtain relative accurate result by considering above features, ASAP-2D (Accident Simulation on Pressure vessel-2 Dimensional) is developed in Shanghai Jiao Tong University to evaluate safety margin of IVR strategy. This paper presents the development of a transient analysis module based on two dimensional spherical coordinate with improved heat conduction model for molten vessel. In additional, POST-CHF heat transfer characteristic is considered and discussed. This module mainly focuses on heat conduction within lower head wall as well as convection occurs at outer surface of RPV. It can be used for the analysis of temperature distribution and, subsequently, the integrity of lower head. Simulation is based on a given heat flux at inner vessel wall provided by other SA codes, like ATHLET-CD. An assessment of IVR capability of large-scale advanced PWR is provided.
2. Mathematic method
The flow chart of ASAP-2D is shown in Fig. 2, where t is time, N is iteration step. Detailed information for calculating heat transfer coefficient (HTC) at outer vessel wall is discussed in Section 2.3.
Fig. 2. Code flow chart.
Fig. 3. Scheme of control volume nodalization.
Control volume scheme are shown in Fig. 3. Nodalization of lower head could be adjusted by user, currently 90 nodes in both r direction and Q direction is applied. Heat conduction between every adjacent two control volume can be obtained according to two-dimensional Fourier equation with finite difference method. Then these combined equations are solved implicitly through Gauss-Seidel iteration.
Equations can be obtained for each node subject to these assumptions:
(1) Research object is lower head with the range from Q = 0°to Q = 90°. And adiabatic boundary condition is applied at Q = 0°and Q = 90°.
(2) For the melted part of the vessel wall, simply modified (Kolev, 2009) specific heat at constant pressure is applied, which means cp is obtained by latent heat of solidification divided by temperature difference between solidus temperature and liquids temperature of vessel.
(3) Radiation heat transfer between vessel wall and other structures in RPV is negligible.
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
2.1. Governing equations
As discussed above, heat conduction equation based on 2-D spherical coordinate is more close to reality. It is more suitable to describe transverse heat conduction effect as well as surface area difference between inner vessel wall and outer vessel wall. Thus RPV heat conduction equation is given by:
dT(r, 0) 1 ^ ( 2 @T(r, 0) r2 dr Vk dr
r2 sin 0
_d_ d0
— IX sin 0
dT(r, 0) d0
dT (r, 0)
dT (r, 0)
= qtn(0), 0 6 0 6 2
= hout(0)-(T(rout, 0)-Tf (0)),
0 6 0 6 p
X dTC^ii = 0,
0 ^ 0; p
Table 2
Heat transfer correlations used in this study.
Correlation
Nu = 0.825 +
[1+(0.492/Pr)a/
Nu = 0.076Ra1/3
ERI (Churchill-Chu) (Churchill and
Chu, 1975) DOE (Churchill-Chu) (Theofanous et al., 1996)
Eckert-type (Eckert and Jackson, 1950) nu = 0.508Pr1/4g0 + Pr) —'[/4Ra,/4 Effective thermal conductivity (He 167 (W/m K) et al., 2005)
where q is vessel wall density in kg/m3, c is specific heat capacity at constant pressure in J/kg K, T is local vessel wall temperature in K, k is local thermal conductivity in W/m K, h is polar angle from bottom center of the lower head in radian, t is time in s.
The boundary conditions for calculation are:
length of molten vessel in h direction is much larger than that in r direction, correlations for natural convection with vertical or inclined cooled surface are used.
It should be noted that Churchill-Chu correlation is suitable for vertical plate and applicability range is any Pr and 0.1 < Ra < 1012, while Eckert-type correlation is suitable for inclined cooled surface. Ra and Pr are defined as follows:
gßL3(T,n - Tmelt)
Pr = -a
where qin (h) is heat flux profile at inner wall in W/m2, rin and rout are inner radius and outer radius of vessel respectively in m, hout (h) is heat transfer coefficient at outer wall, which depends on local thermal physical condition in W/m2 K, Tf (h) is local fluid temperature in ERVC channel in K.
where L is the width of molten vessel in m, Tin is wall temperature at inner surface in K, Tmelt is melting temperature in K.
For wall temperature between solidus temperature and liquids temperature, linear interpolation of solid vessel conductivity and molten vessel conductivity is applied.
2.3. Heat transfer logic at the outer surface
2.2. Thermal conductivity of the wall
It is known that thermal conductivity would affect temperature field of lower head. In this study, when wall temperature is lower than solidus temperature, thermal conductivity is obtained according to Theofanous et al. (1996) shown in Fig. 4. However there is little research available for thermal conductivity when wall temperature higher than liquids temperature. Thus it is assumed that molten vessel stay in place and local natural convection could occur. This study attempts to capture some effects of convection within the conduction model. Both natural convection correlation and effective thermal conductivity method (Dinh et al., 1996; He et al., 2005) is studied. Empirical natural convection heat transfer correlations and effective thermal conductivity are listed in Table 2. Since the
o n TD £= O O
Temperature (K)
Fig. 4. Thermal conductivity for temperature lower than solidus temperature.
Since external vessel heat transfer models could be in different kind of region, which depends on local thermal hydraulic behavior. Within this code, heat transfer mode at outer surface wall consists of liquid phase convection, subcooled nucleate boiling, saturated nucleate boiling, transition boiling and film boiling. Heat transfer correlation logic for above mentioned correlation is shown in Fig. 5.
Tonb is wall temperature corresponding to onset nucleate boiling obtained by Bergles and Rohsenow correlation (1963). The equation is expressed in Eq. (7):
Tonb = Tsat + 0.556
qout(0)
1082p1
0.463p00234
where Tsat is saturated fluid temperature in K, qout (h) is heat flux at outer surface wall in W/m2, p is pressure in bar. Hf is ERVC fluid enthalpy in J/kg, Hsat is saturated liquid phase enthalpy in J/kg, qchf is critical heat flux obtained by correlations in W/m2, Tchf is wall temperature corresponding to qchf in K.
Tchf is reference wall temperature when heat flux reaches CHF, in
Tmfb is minimum film boiling temperature obtained by Groeneveld and Stewart's correlation (1982), which is also used in TRACE, CATHARE and ATHLET. The correlation is:
For P 6 9 x 106 (Pa):
Tmfb = 557.85 + 0.0441
- max 0.0
104xq,
2.82 + 0.00122 x 10-3P
For P >9 x 106 (Pa):
TMFB = Tsat +(TMFB,9 Mpa — Tsat-^-7—6
Pc — 9x10
r=rout
10—33P- 3.72 x 10—12P2
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
Fig. 5. Heat transfer mode logic.
Table 3
Heat transfer model.
Regimes Model
Liquid phase convection Dittus-Boelter
Subcooled nucleate boiling Modified Chen
Saturated nucleate boiling Chen
Transition boiling Interpolation between nucleate boiling and
film boiling
Film boiling Berenson
Nu = 0.023 • Re0 8 ■ Pr'
divided into latent heat for vapor formation, convection caused by bubble agitation and single-phase heat transfer, thus:
Qout = Qe + Qa + Qs,
where p is pressure in Pa, xqm is steam mass quality, pc is critical pressure in Pa, TMFBr 9Mpa is TMFB at 9 MPa in K, Tsat is saturated fluid temperature at 9 MPa in K.
2.4. Heat transfer correlation at the outer surface
Wall-to-fluid heat transfer mode is selected from representative computer codes that are most widely used, such as ATHLET, RELAP5, MELCOR, TRAC-M, and COBRA. Heat transfer models are summarized in Table 3 (The SCDAP/RELAP5-3D Code Development Team, 2003; Austregesilo, 2012; Collier and Thome, 1994; Choi et al., 2009). The correlations listed in Table 3 will be discussed in the following subsection.
2.4.1. Liquid phase convection
For wall temperature lower than Tonb, liquid phase convection mode is used. In current study, the inlet mass flow rate in ERVC channel is time dependent. Since Reynolds number in ERVC case is normally larger than 2300, the well-known Dittus-Boelter correlation is selected for liquid phase convection, which is given by:
Then ''pumping factor'' is defined by the ratio of heat flux due to bubble agitation and that causing vapor formation (Hari and Hassan, 2002) as:
e = (13)
After the model is compared to available data, it is assumed that parameter e is piecewise constant:
T3.2 •
e = < 1.3
prcp qv -(hvs-hls)
0.1 Mpa < P < 0.95 Mpa 0.95 Mpa < P < 5 Mpa 5 MPa < P
where ql and pv are density of liquid and vapor, respectively, kg/m3, cp is specific heat capacity, J/kg K, hts and hvs are saturation enthalpy of liquid phase and vapor phase, respectively, J/kg. Then local steam mass quality comes from:
Xqm (h) =
j (Qout(h)~ Qsp (h)) • Aout(0) Gj •(hvs - hls)• (1 + e)
with j is control volume number, xqm is steam mass quality in control volume j, Nonb is the number of first nucleate boiling control volume, Aout is heat surface area, m2, Gj is ERVC mass flow rate, kg/s, qsp is given by:
Qsp(h) =
J 0.0 j
\Qout W-hh-^j
: < Nonb P Nonb
where honb is enthalpy corresponding to onset nucleate boiling temperature.
2.4.2. Subcooled nucleate boiling
If wall temperature is higher than Tonb and fluid remain in sub-cooled boiling state, modified Chen correlation is used. It is assumed that total heat flux is composed of nucleate boiling part and single-phase convective part:
qout = hNB(T - TSat)+ hc(T - Tfluid) (11)
where hc and hNB are heat transfer contribution due to convection and nucleate boiling, respectively. Detailed information is provided by John G. Collier and Thome (1994). In order to calculate hNB and hc, steam mass quality is obtained according to Bowring model (Collier and Thome, 1994), in which heat removed from vessel wall is
2.4.3. Saturated nucleate boiling
Chen correlation (1966) is adopted in this code when wall temperature is higher than Tonb, and ERVC fluid reaches saturated state. Total heat flux transfer coefficient at outer surface wall can be expressed as follows:
hchen = hmac • F + hmic • S (17)
where hmac and hmic are the macroscopic and microscopic heat transfer coefficient respectively. Specifically, hmac is obtained by convec-tive heat transfer correlation and hmin is calculated from Forster and Zuber correlation (1955). F and S stand for Reynolds number factor and suppression factor, respectively.
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
5 J 1.5 1.0 0.5
position (degree)
С 10-
ф о iiz Ф
J» 10^
(Л £= ГО
Is 101
£ 400
RELAP5
linear(SCDAP)
ATHLET(COSINE)
MELCOR(LOG)
CATHARE
500 600 700 Temperature (K)
Fig. 6. Comparison with different CHF correlation prediction.
Fig. 7. Transition boiling models.
2.4.4. Critical heat flux
Compared to ERVC case, most of CHF correlations are validated in relative high pressure. Within correlations database provided by Todreas and Kazimi (2010), only Biasi and Bowring CHF correlation databases cover 0.2 MPa condition. However CHF predicted by these two correlations are significant different from that of Theofanous et al. (1996) correlation. As is shown in Fig. 6, CHF obtained by Bowring correlation is overestimated, while CHF obtained by Biasi correlation is much lower than that of Theofanous. Thus correlation from Ttheofanous is applied which is given as:
qCHF(h)=490 + 30.20 - 8.88 • 10-1h2 + 1.35 • 10-2h3 - 6.65
• 10-V
where в is in degree and qchf is in kW/m2
"2 ф œ E ф
2.35 2.40
2.45 2.50 thickness (m)
2.55 2.60
Fig. 8. Verification result.
2.4.5. Transition boiling
When wall temperature exceeds Tchf or heat flux is higher than qCHF, POST-CHF correlation is applied. If local wall temperature is lower than Tmfb, transition boiling correlation is adopted. Since there is no widely accepted transition boiling model, most existing safety analysis codes apply simplified correlations which consider contribution by nucleate boiling and by film boiling. Different weighting function is summarized in Table 4.
Calculation results by these candidate transition boiling models are shown in Fig. 7. In this case, heat flux predicted by ATHLET and SCDAP show highest value among these codes. Heat flux obtained by CATHARE and TRACE is similar but they produce a little smaller heat flux than ATHLET. Transition heat flux predicted by MELCOR
Table 4
Summary of transition boiling model.
Transition boiling model
ATHLET MELCOR
TRACE (Choi et al., 2009)
CATHARE (Choi et al., 2009)
Cosine shape interpolation (Austregesilo, 2012) Logarithmic interpolation (Gauntt, 2005)
qTB = SqcHF + (1 -d = p1 - rti!--TM" '
- d)qFB
RELAP5/MOD3 SCDAP/RELAP
\Tchf -tmfb
Tmfb : Groeneveld-Stewart's correlation Jones's model qTB = SqCHF + (1 - S)qFB
d = ( tw-tmfb\ 2 \tchf -tmfb )
TMFB: Groeneveld-Stewart's correlation
Chen-Sundaram-Ozkaynak (The RELAP5 Code Development Team, 1995)
Linear interpolation (Ahn and Kim, 2003)
0.7 0.6 g 0.5
Il 0.4
ë 0.3 го
ë 0.2 0.1
-■— DOE
ASAP-2D-1
ASAP-2D-2
ASAP-1D-2
ASAP-1D-1
angle (o)
Fig. 9. Comparison of ASAP outer vessel heat flux with DOE.
decreases faster than the three models mentioned above. Chen-S undaram-Ozkaynak model in RELAP5/MOD3 produces smallest heat flux among these codes. Since Chen-Sundaram-Ozkaynak model deviates far from other models, it is not applied in current version of code. Simulation results show that, compared with other heat transfer region, transition heat flux obtained by different codes varies significantly.
2.4.6. Film boiling
Film boiling region is selected if wall temperature is higher than Tmfb. In current case, only Berenson model (Collier and Thome,
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
Table 5
Key input parameter.
Parameter Value
External fluid temperature (K) 400
Thermal conductivity (W/m K) 32
Lower head Inner radius (m) 2
Melting temperature (K) 1600
Thickness (m) 0.15
Volume (m3) 10.85
Density (kg/m3) 8191
Molten ceramic pool Height of pool (m) 1.52
Pool angel (°) 76.14
Heat generation (MW/m3) 1.3
Density (kg/m3) 6899.2
Metallic pool Height of pool (m) 0.9273
Heat generation (MW/m3) 0
angle ( )
Fig. 11. Comparison of output heat fluxes and qchf.
ä 10 H
m <0 £=
- DOE ASAP-2D-1 ASAP-2D-2 ASAP-1D-2 ASAP-1D-1
angle (o)
Fig. 10. Predictions for vessel wall thickness.
angle ( )
Fig. 12. Comparison of CHF ratios.
Table 6
Input parameter for MOPOL.
Parameter (unit) Value
1000 MWe class 1700 MWe class
Mass of stainless steel (t) 28.4 Proportion of melting UO2 (%) 86.6 Decay heat (MW) 27.1 Proportion of oxidized Zr (%) 39.5 Lower head radius (m) 2.0 Lower head thickness (m) 0.15 Gap between vessel and insulation (m) 0.15 68.8 86.6 46.1 39.5 2.523 0.192 0.15
Table 7 Boundary conditions and heat transfer models.
Item Content
Melt pool heat transfer correlations External vessel mass flow rate External vessel inlet water temperature External vessel inlet water pressure Effective conductivity for molten vessel Interpolation method for transition boiling Mini-ACOPO 452 (kg/s) 323 (K) 0.18 (MPa) DOE (Churchill-Chu) Logarithmic interpolation
hFB = 0.425
g -{Pl - Pv)-Pv • 4 • h'lv
r{pl-pv )
2.0x104 1.5x104 1.0x104 5.0x103 0.0
N . —»— heat transfer coefficient
\ —o— vessel thickness
onset of \ ^^^ '
nucleate
boiling \ "
0.20 0.15
<f) a>
0.05 S
angle (O)
Fig. 13. Heat transfer coefficient and vessel thickness.
where g is gravity acceleration, in m/s2, kv is thermal conductivity for vapor, in W/m K, iv is dynamic viscosity for vapor, in Pa s, DT is super heat, in K, r is surface tension, in N/m, h'lv is an effective latent heat of vaporization allowing for the effect of superheat:
1994) is applied, which incorporate hydrodynamic wave instability theory of Taylor into the Bromley model:
h'lv = hlv •
Cp • DT
3. Verification
In order to evaluate feasibility of the new developed module, simulation results are compared to one-dimension analytical solution with uniform heat flux distribution.
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
1800 1900
2,5x106 2,0x106 1,5x106 1,0x106 5,0x105 0,0
—■— Qchf 1 1
acopo qin acopo qout
k-m qout Ref. qin Ref. qout
Fig. 14. Temperature field at 5127s.
angle ( )
Fig. 17. Heat flux distribution.
0 1000 2000 3000 4000 5000 6000 Time (s)
Fig. 15. Transient temperature distribution at inner vessel wall.
angle (
Fig. 18. Vessel thickness distribution.
2.0x104
■ OUT 1
• OUT 10
A OUT 20
T OUT 30
< OUT 40
OUT 50
♦ OUT 60
• OUT 70
• OUT 80
» OUT 90
0 2000 4000 6000
Time (s)
Fig. 16. Transient temperature distribution at outer vessel wall.
1.5x104 -
1.0x10-
£= CD O Ú= CD O O
T« 5.0x103-| £= TO
angle ( )
Fig. 19. External vessel heat transfer coefficient distribution.
Since it is difficult to obtain analytical solution based on two-dimensional spherical coordinate with practical boundary condition, one-dimensional analytical solution with uniform heat flux input is applied in this section. Consider a hollow sphere with inner and outer radii rin and rout, respectively. It is assumed that Neumann (prescribed heat flux as boundary condition) and Robin (prescribed heat transfer coefficient and fluid temperature) boundary conditions are applied at inner surface and outer surface, respectively. Both thermal conductivity and heat transfer coefficient is assumed to be constant in the section. At steady state, heat conduction in spherical coordinates becomes:
dr (r dr) 0
The general solution of which is:
T = — + c2 r2
Since Nuemann boundary condition is applied at inner surface of vessel, for r = rin has:
Q(rin) = -k dr
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
£ 30-1
Film boiling
occurs at 75
Film boiling occurs at 2942s
5000 10000 Time (s)
Fig. 20. POST-CHF zone spread (qin is obtained according to ACOPO correlation).
0 10 20 30 40 50 60 70 80 90 angle (O)
Fig. 23. External heat transfer coefficient.
2.5x106 _ 2.0x106-
$ 1.5x106-x
S 1.0x10-ro
c 5.0x105-0.0
—■— input heat flux 1
200 kg/s
Qchf -
700 kg/s
1000 kg/s
angle ( )
Fig. 21. External cooling mass flow rate.
2.5 2.0 1.5 1.0 0.5 0.0
- input heat flux Qout(Ref.) Qchf
Qout(T=60 oC) Qout(T=70 oC) Qout(T=80 °C)
angle ( )
Fig. 24. Heat flux distribution.
v> m cu £=
yielding:
r2q{rin)
-«— 200 kg/s
—— Ref.
700 kg/s -
—1000 kg/s
angle ( )
Fig. 22. Vessel thickness.
Robin boundary condition leads to:
r¡nq(rm)
r20uth
+ tf -
so that the unknown temperature at r = rx is:
T(r, rlq{rin) r2nq{rin) t rinq{rin)
T{rx)=~ir— + + tf
r2cuth
E 0.15-
v> m cu £=
angle (O)
Fig. 25. Vessel thickness distribution.
Fig. 8 compares the analytical result and numerical result, where h = 10,000 W/m2 is heat transfer coefficient at outer surface of vessel, k = 32 W/m K is thermal conductivity of vessel, tf =373 K is external vessel fluid temperature, q (rin) = 1 x 106 W/m2.
Good agreement between the prediction and analytical result indicates that new developed module is feasible to simulate one-dimensional temperature distribution with uniform heat flux input. The module development is free of obvious errors and mistakes.
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
£= CD O
£ m £= TO
—«— Ref.
—— T=333K
T=343K
T=353K
angle ( )
Fig. 26. External vessel heat transfer coefficient distribution.
2,5x106-|
2,0x10-
§ 1,5x10-
^ 1,0x106 -| to
^ 5,0x105 -
—■— input heat flux ' i '
cosine(ATHLET-CD •
linear(SCDAP) -
log(MELCOR)
—CATHARE
. 'jf -
angle ( )
Fig. 27. Heat flux distribution.
4. Application
In this section, this module is applied to final bounding state (FIBS) benchmark result conducted by DOE and large scale advanced PWR.
4.1. Applied to benchmark calculation
First this module is used to benchmark calculations with input parameter in NUREG/CR-6489, in which heat transfer behavior of melt pool and vessel wall of AP600 is studied (Esmaili and Khatib-Rahbar, 2004). In this benchmark, heat flux imposed on lower head is predicted by different correlations for melt pool. Also thickness of crust and vessel wall is analyzed. Since only one-dimensional heat conduction in the vessel wall is considered in DOE benchmark, DOE output heat flux distribution is the same with the input. In order to study two dimensional effect brought by this typical non-uniform heat flux, heat flux derived from DOE benchmark is applied as boundary condition for ASAP-2D.
Key input parameters are summarized in Table 5. In this section, four different models of ASAP-2D code are adopted for comparison. Then vessel thickness and heat flux distribution at outer vessel is compared to benchmark.
Figs. 9 and 10 show comparison of benchmark results and simulation results obtained by ASAP:
ASAP-2D-1: two-dimensional heat conduction and constant
thermal conductivity.
ASAP-2D-2: two-dimensional heat conduction and thermal conductivity vary with temperature (shown in Fig. 4). ASAP-1D-1: one-dimensional heat conduction and constant thermal conductivity.
ASAP-1D-2: one-dimensional heat conduction and thermal conductivity vary with temperature (shown in Fig. 4).
It can be seen that results obtained by ASAP have good qualitative agreement with benchmark results. But heat flux is notably lower than benchmark due to surface area difference and/or two-dimensional heat conduction effect. As is illustrated in Fig. 9, both ASAP-1D-1 and ASAP-1D-2's predictions are almost parallel to DOE's result but smaller than it. This is due to outer surface area is larger than inner surface. ASAP-2D-2's prediction in Fig. 9 indicates that heat flux distribution is flatten, especially at location adjacent to the metallic layer (larger than ~73°), when two-dimensional heat conduction is considered.
Fig. 10 compares vessel wall thickness predicted by DOE and ASAP. Except for ASAP-2D-2, all the other three results obtained by ASAP are slightly larger than DOE. This is because output heat flux is lower than DOE. ASAP-2D-2 prediction is notable larger than DOE due to thermal conductivity varies with temperature. As is shown in Fig. 4, thermal conductivity is larger than DOE when wall temperature is lower than about 900 K. Since wall temperature is relatively low at location adjacent to ceramic pool (smaller than ~72°), transverse heat conduction is enhanced. Thus, ''focusing effect'' is further alleviated, which leads to larger solid vessel thickness.
It can be seen from Figs. 9 and 10 that ASAP shows reasonable applicability to simulate the thermal behavior of pressure vessel during severe accident. After taking 2-D heat conduction into consideration, maximum heat flux could decrease 23% and minimum vessel thickness could increase 60%. Results indicate that 2-D heat conduction effect plays an important role in analyzing temperature field of lower head. Thus this effect should be considered. Since the prediction is more accord with reality, two-dimensional heat conduction model and thermal conductivity varying with temperature model is used in the following work.
4.2. Applied to large-scale advanced pressurized water reactor
One of most challenging problems for large-scale advanced pressurized water reactor is the feasibility of IVR strategy. Since simulation model provided in this paper is more close to reality, the relatively accurate temperature field could help in identifying the integrity of lower head. It is possible that safety margin could be increased to some extent.
Both RPV size and heat flux non-uniformity would be enlarged in large-scale advanced pressurized water reactor. Since two-dimensional heat conduction effect would play a more important role when heat flux non-uniformity or lower head radius becomes larger, it is important to analysis integrity of lower head. Thus preliminary thermal analysis on lower head of a 1700 MWe class advanced PWR during postulated severe accident is performed.
4.2.1. Boundary condition
Since detailed design parameter is not complete yet, it is hard to set up an accurate 1700 MWe class reactor model. As a boundary condition, heat flux load is obtained by a lumped parameter code MOPOL, which is developed in Shanghai Jiao Tong University for simulating melt pool behavior at final bonding state. Key input parameters for MOPOL are listed in Table 6, where the value of 1000 MWe class reactor is derived from a bounding severe accident analysis based on MAAP4 (Bao, 2013). Due to 1700 MWe class advanced PWR adopts similar design with that of AP1000, it is
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
assumed that same fuel assembly structure is applied. Thus decay heat and mass of stainless steel is scaled from 1000 MW class reactor to 1700 MW class reactor. Although thermal power of 1700 MWe class reactor is larger than that of AP1000, CHF correlation for AP1000 is applied this part for conservative reason. The CHF correlation is given as (Esmaili and Khatib-Rahbar, 2005):
qCHF(h) = 1.44 x(490 + 30.20 - 8.88 • 10-1h2 + 1.35 • 10-2
• 03 - 6.65 • 10-5h4) (27)
where 0 is in degree and qchf is in kW/m2.
Currently, two-layer molten pool configuration is selected for melt pool. Simulation is terminated when energy difference between inner surface and outer surface is less than 0.5% or total calculation time reaches 15,000 s. Boundary conditions and heat transfer models applied in reference case are summarized in Table 7. (Gauntt, 2005; Choi et al., 2009; Bao, 2013; Esmaili and Khatib-Rahbar, 2005) Some of these items are chosen to perform sensitivity analysis later.
4.2.2. Reference case calculation results
The simulation results of in-vessel retention are illustrated in this section. Calculation results are compared with MOPOL results. Since one-dimensional heat conduction model in vessel wall is applied in MOPOL, output heat flux distribution is the same with input one. It can be seen from Figs. 11 and 12 that ASAP output heat flux shows similar trend with MOPOL output heat flux. Meanwhile, heat flux predicted by ASAP is significantly lower than that of MOPOL due to two-dimensional heat conduction effect and surface area difference. In MOPOL prediction, CHF ratio (local output heat flux divided by local CHF) reaches maximum (0.91) at the location near 69.8°, where is the interface of ceramic pool and metallic layer. However, maximum CHF ratio predicted by ASAP (0.78) is about 15% lower than MOPOL prediction.
Heat transfer coefficient at external vessel wall and vessel thickness is shown in Fig. 13. With increasing polar angle, heat transfer coefficient decreases until reach onset of nucleate boiling point (~19°). This is because before onset of nucleate boiling point, heat transfer mode is single phase convection. Increasing cross section area of external reactor cooling channel leads to decreasing flow velocity. Then subcooled boiling occurs and heat transfer coefficient begins to increase due to increasing convective heat transfer effect brought by overheated vessel wall. The vessel wall is molten to some extent for the position between 17°and 90°. While between 60° and 90°, the thickness of solid vessel wall is no more than 3 cm. This is because input heat flux in this area is relatively larger than others, which leads to larger temperature gradient and smaller wall thickness. Temperature field at final time point (5127 s) is illustrated in Fig. 14. The black line in Fig. 14 shows melting temperature of the lower head. At the end of simulation, nodes 20-90 at inner vessel reach melting temperature. However, temperature of all the nodes at outer vessel wall is lower than 450 K. This is because external heat transfer mode in this zone is subcooled nucleate boiling, which could provide sufficient cooling.
Transient temperature distribution at inner vessel wall and outer vessel wall are shown in Figs. 15 and 16. The position node 1 is at the bottom center of lower head and node 90 is at top of metallic layer. It can be seen that temperature for node 50 to node 90 reach melting temperature soon after simulation begins. This is due to heat flux imposed on these nodes is relatively higher than others.
It can be seen that output heat flux predicted by ASAP is relatively smaller than MOPOL prediction. Thus higher margin of safety could be obtained. POST-CHF heat transfer does not occur due to all the nodes at outer surface are lower than 450 K. However the thick-
should be done to examine whether 3 cm is thick enough to sustain the weight of lower head. Otherwise, improvement should be done to the design of large-scale advanced pressurized water reactor, such as add sacrificial material to increase the thickness of metallic layer aiming to mitigate ''focusing effect''.
5. Sensitivity analysis
A set of sensitivity analysis simulations are performed in this section to study the influence of boundary conditions and heat transfer correlations. For boundary conditions, different melt pool heat transfer correlations, external reactor vessel cooling mass flow rate and external vessel inlet coolant temperature are taken into consideration. Different transition boiling correlations and effective conductivity correlations for molten vessel are also chosen to perform sensitive analysis. Boundary conditions and heat transfer models used in reference case are listed in Table 7.
5.1. Me't poo' heat transfer correlations
In this section, two sets of heat transfer correlations for melt pool are chosen to perform sensitivity analysis. One set is K-M, which means that heat transfer correlations for top and bottom surface of ceramic pool are Kulacki-Emara (1975) and Mayinger et al. (1976). The other set is ACOPO (Esmaili and Khatib-Rahbar, 2005), which means correlations for top and bottom surface of ceramic pool are derived from ACOPO experiment. The correlation applied in reference case is mini-acopo correlation. In the results for AP600, heat flux predicted by K-M and mini-acopo shows good agreement with each other (Esmaili and Khatib-Rahbar, 2004). Since in this large-scale advanced PWR both radius of RPV and thermal power is increased, Rayleigh number (based on internal heat generation) are of the order 1017. Thus in K-M and mini-acopo correlations, different ratio of heat which is transferred to top and lower boundary of melt pool could be found. Then this would lead to a notable difference in heat flux distribution.
Results are shown in Figs. 17-19. qin show heat flux imposed on inner vessel wall and qout shows heat flux distribution at outer vessel wall. In Fig. 17, it can be seen that qout could decrease to 0.25 MW/m2 when ACOPO or K-M correlations are applied. This is because film boiling correlation is chosen when wall temperature is larger than TmSb. As a result, heat transfer coefficient decrease to about 150 W/m2 K at location from 45°to 90° (in Fig. 19). Thus heat conduction in h direction is enhanced which leads to decreasing qout in this zone. Vessel thickness is illustrated in Fig. 18, for the location larger than 45°, vessel wall is completely molten due to large input heat flux and lack of effective cooling. Besides, qin is larger than qchf only in part of nodes (between 70°and 80°) while POST-CHF heat transfer occurs at the location between 45° and 90°. This is because heat flux in r direction (Fig. 1) would decrease rapidly when POST-CHF heat transfer region occurs in any node. Thus transverse heat conduction in h direction is enhanced, which leads to increasing temperature of nodes adjacent to it. POST-CHF zone would spread until temperature of following nodes below Tchf. In Fig. 20, film boiling region, shown in shaded area, first occurs in node 75 at 2942 s. Then POST-CHF zone continually spread due to transverse heat conduction in h direction.
Besides, CHF ratio would be still less than unity when POST-CHF heat transfer occurs and leads to the failure of lower head. This is because decreasing external heat transfer could enhance transverse heat conduction, leading to decreasing in heat flux in r direction in POST-CHF region. Thus it should be noted that when taking two-dimensional heat conduction into consideration, lower head
ness of lower head is no more than 3 cm at 60-90°. Further study may fail even when CHF ratio is less than unity.
CATHARE trace
log(MELCOR) linear(SCDAP) cosine(ATHLET-CD)
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
angle ( )
Fig. 28. Vessel thickness.
J=. 0.1505
Ü 0.10-1 o
$ 0.05-1
—■— constant X=167 w/mk
—•— Ref.
eckert
constant X=32 w/mk
angle ( )
Fig. 31. Vessel thickness.
1051—p
c 104^
- cosine(ATHLET-CD) " linear(SCDAP) log(MELCOR) trace
CATHARE
angle ( )
Fig. 29. Heat transfer coefficient.
£= CD O à= CD O O
m £= to
—■— constant X=32 w/mk 1 '
■ eri
eckert
. constant X=167w/mk
angle ( )
Fig. 32. External heat transfer coefficient.
2.5x106-_ 2.0x106-
§ 1.5x106 -
t; 1.0x106 -
1 5.0x105-
0.00 30 60 90 angle (O)
Fig. 30. Heat flux distribution.
5.2. External cooling mass flow rate
Natural circulation flow behavior is also an important factor affecting heat removal capacity in IVR strategy. Since constant mass flow rate for ERVC channel is applied as boundary condition in current version, sensitive analysis on external cooling mass flow rate should be considered. In this part, three different value of mass flow rate are chosen to perform sensitivity analysis. Results are shown in Figs. 21-23. It can be seen from Fig. 21 that output heat flux is
similar to each other. The only difference is that, at the location between 0°and 42°, heat flux is slightly larger than other prediction using 1000 kg/s. From Fig. 22 it can also be found that vessel thickness begins to decrease at 13°, which is apparently different from others. This is because heat transfer coefficient (Fig. 23) is larger than others at the location between 0°and 42°.
5.3. External vessel inlet coolant temperature
In this section, sensitivity to external vessel inlet coolant temperature is studied. Coolant inlet temperature is set to 333, 343 and 353, respectively. In order to focus on this effect, it is assumed that mass flow rate is the same to reference case. Results are illustrated in Figs. 24-26. It can be seen from Fig. 24 that external heat flux is almost the same except for the node smaller than 45. As is shown in Fig. 25, vessel thickness begins to decrease at node 15 when inlet coolant temperature is set to 343 and 353 K. This is because output heat flux between node 15 and 45 is relatively larger in these two cases. In Fig. 26, we can see that with increasing coolant inlet temperature, external heat transfer coefficient is also increased. This is because nucleate boiling effect increases with increasing inlet coolant temperature.
5.4. Transition boiling correlation
In order to examine sensitivity to transition boiling correlation, ACOPO correlations for melt pool are used to calculate input heat flux. Five different transition boiling correlations from typical
Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
system codes (shown in Table 4) are chosen to perform sensitive analysis. Figs. 27-29 show the distributions of heat flux, vessel thickness and external heat transfer coefficient. These results can be categorized into two groups according to whether vessel wall is molten through or not. Results in the same group are similar to each other while obviously different from the cases from the other group. It can be seen that TRACE or MELCOR transition boiling heat transfer correlation would lead to the failure of lower head. At the nodes larger than 40, transition boiling turn into film boiling rapidly due to fast increase in wall temperature. However reactor vessel does not fail in other three cases. This is because transition boiling heat transfer correlation among these codes is different from each other. As is shown in Fig. 7, MELCOR and TRACE produce lowest transition boiling heat transfer coefficient, which means wall temperature would increase faster than others. Then transverse heat conduction leads to rapid temperature increase for adjacent nodes. Thus POST-CHF zone is easier to spread. On the contrary, since wall temperature in this group is usually less than 460 K, transition boiling heat transfer coefficient predicted by ATHLET, SCDAP/RELAP and CATHARE is relatively high, which mitigates wall temperature rise in transition boiling region. Thus reactor vessel is less likely to fail.
5.5. Effective conductivity correlation
Since thermal conductivity affects temperature field of lower head, different heat convection correlations and effective thermal conductivity are applied when local wall temperature exceeds liquids temperature (shown in Table 2). In order to study the difference among these correlations and effective thermal conductivity, sensitivity analysis is performed in this section. In reference case, DOE (Churchill-Chu) correlation is applied. And results are also compared with the predictions on base of constant thermal conductivity (32 W/m K is derived from NUREG/CR-6489). Results are shown in Figs. 30-32.
It can be seen that DOE, ERI and Eckert-type correlation produce a very similar results. But these results are different from constant thermal conductivity model (32 or 167). This is because three correlations produce a higher heat transfer coefficient at the location between 15 and 45, which leads to higher heat flux and lower vessel thickness.
6. Conclusion
A two dimensional transient analysis module, namely ASAP-2D, for integrity of reactor lower head during severe accident is developed. The module is verified by one dimensional analytical result. Then it is applied to DOE benchmark and large scale advanced pressurized water reactor. Conclusions are listed as follows:
(1) One-dimensional temperature distribution with uniform heat flux input predicted by ASAP shows good agreement with analytical result. This indicates that ASAP is suitable to simulate heat transfer characteristic of lower head during severe accident.
(2) Compared to DOE benchmark, ASAP produces similar heat flux distribution and vessel thickness distribution. Meanwhile after considering two-dimensional effect, maximum output heat flux could decrease 23%. Thus it is recommended to consider 2-D effect when IVR analysis is performed.
(3) Sensitivity analysis is carried out for boundary conditions and new developed models, including heat transfer
correlations for melt pool, mass flow rate in ERVC, transition boiling heat transfer correlation, effective conductivity correlation for molten vessel and ERVC inlet water temperature. Results indicate that for large-scale advanced PWR outer heat flux distribution is more sensitive to heat transfer correlations for melt pool (although the results in AP600 are not very sensitive to these correlations) and transition boiling correlation than the others.
(4) CHF ratio should be carefully examined when it is used as criteria for lower head failure. When taking two-dimensional heat conduction into consideration, in some cases, lower head may fail even when CHF ratio is less than unity.
Acknowledgement
The authors would like to thank GRS for providing the ATHLET-CD code.
References
Ahn, K.I., Kim, D.H., 2003. A state-of-the-art review of the reactor lower head models employed in three representative US severe accident codes. Prog. Nucl. Energy 42 (3), 361-382.
Austregesilo, H., 2012. ATHLET-CD mod 2.2 - cycle C User's manual.
Bao, H., 2013. Study on Heat Transfer Models of Molten Pool in the Lower Head during a Severe Accident. Master thesis, Shanghai Jiao Tong University, P.R. China.
Bergles, A.E., Rohsenow, W.M., 1964. The determination of forced-convection surface-boiling heat transfer. J. Heat Transf. 86 (3), 365-372.
Chen, J.C., 1966. Correlation for boiling heat transfer to saturated fluids in convective flow. Ind. Eng. Chem. Process Des. Dev. 5 (3), 322-329.
Choi, K.Y., Yun, B.J., Park, H.S., Kim, H.D., Kim, Y.S., Lee, K.Y., Kim, K.D., 2009. Development of a wall-to-fluid heat transfer package for the SPACE code. Nucl. Eng. Technol. 41 (9), 1143-1156.
Churchill, S.W., Chu, H.H., 1975. Correlating equations for laminar and turbulent free convection from a vertical plate. Int. J. Heat Mass Transf. 18 (11), 1323-1329.
Collier, J.G., Thome, J.R., 1994. Convective Boiling and Condensation. Oxford University Press.
Dinh, T.N., Bui, V.A., Nourgaliev, R.R., Sehgal, B.R. 1996. Crust dynamics under PWR in-vessel melt retention conditions. In: ANS Proc. of 1996 National Heat Transfer Conference, Texas, 368-375.
Eckert, E.R., Jackson, T.W., 1950. Analysis of Turbulent Free-Convection Boundary Layer on Flat Plate (No. NACA-TN-2207). National Aeronautics and Space Administration, Washington DC.
Esmaili, H., Khatib-Rahbar, M., 2004. Analysis of In-Vessel Retention and Ex-Vessel Fuel Coolant Interaction for AP1000, NUREG/CR-6849.
Esmaili, H., Khatib-Rahbar, M., 2005. Analysis of likelihood of lower head failure and ex-vessel fuel coolant interaction energetics for AP1000. Nucl. Eng. Des. 235 (15), 1583-1605.
Forster, H.K., Zuber, N., 1955. Dynamics of vapor bubbles and boiling heat transfer. AIChEJ. 1 (4), 531-535.
Gauntt, R.O., 2005. MELCOR Computer Code Manual Version 1.8.6.
Groeneveld, D.C., Stewart, J.C., 1982. The minimum film boiling temperature for water during film boiling collapse. In: Proc. 7th Int. Heat Transfer Conf., Munchen, 1982.
Hari, S., Hassan, Y.A., 2002. Improvement of the subcooled boiling model for low-pressure conditions in thermal-hydraulic codes. Nucl. Eng. Des. 216 (1), 139-152.
He, X., Elmer, J.W., DebRoy, T., 2005. Heat transfer and fluid flow in laser microwelding. J. Appl. Phys. 97 (8), 084909.
Kolev, N.I., 2009. Multiphase Flow Dynamics 4: Nuclear Thermal Hydraulics. Springer.
Kulacki, F.A., Emara, A.A., 1975. High Rayleigh Number Convection in Enclosed Fluid Layers with Internal Heat Sources. USNRC, NUREG-75/065.
MAAP4, 1994. MAAP4: Modular Accident Analysis Program for LWR Plants, Code Manual Vols. 1-4. Prepared by Fauske & Associates Inc, Burr Ridge, IL, USA for the EPRI, Palo Alto, CA, USA.
Mayinger, F., et al., 1976. Examination of Thermohydraulic Processes and Heat Transfer in a Core Melt. BMFT RS 48/1, Institute für Verfanhrenstechnic der TU, Hanover, Germany.
Rempe, J.L., Knudson, D.L., Allison, C.M., Thinnes, G.L., 1997. Potential for AP600 In-Vessel Retention through Ex-Vessel Flooding. Lockheed Idaho Technologies Co., Idaho National Engineering and Environmental Lab., Idaho Falls, ID (United
14 Z. Cao et al./Annals of Nuclear Energy xxx (2015) xxx-xxx
States). Funding organisation: USDOE Assistant Secretary for Nuclear Energy, Washington, DC (United States). Schulz, T.L., 2006. Westinghouse AP1000 advanced passive plant. Nucl. Eng. Des. 236 (14), 1547-1557.
Sehgal, B.R., Theerthan, A., Giri, A., Karbojian, 2003. Assessment of reactor vessel
integrity (ARVI). Nucl. Eng. Des. 221 (1), 23-53. The RELAP5 Code Development Team, 1995. RELAP5/Mod3 Code Manual. Idaho National Engineering Lab., Idaho Falls, ID, USA.
The SCDAPRELAPS-3D Code Development Team, 2003. The SCDAPRELAPS-3D Code
Manuals, INEEL, Idaho Falls, ID, USA. Theofanous, T.G., Liu, C., Additon, S., 1996. In-Vessel Coolability and Retention of a
Core Melt. D0E/lD-10460 Vol. 1. Todreas, N.E., Kazimi, M.S., 2010. Nuclear Systems. CRC Press. Wittneben, B.B., 2012. The impact of the Fukushima nuclear accident on European energy policy. Environ. Sci. Pol. 15 (1), 1-3.
