The Effective Convectivity Model for Simulation of Molten Metal Layer Heat Transfer in a Boiling Water Reactor Lower Head

Tran, Chi-Thanh; Kudinov, Pavel

doi:https://doi.org/10.1155/2013/231501

Science and Technology of Nuclear Installations

On this page

Abstract Introduction References Copyright Related Articles

Research Article | Open Access

Volume 2013 | Article ID 231501 | https://doi.org/10.1155/2013/231501

The Effective Convectivity Model for Simulation of Molten Metal Layer Heat Transfer in a Boiling Water Reactor Lower Head

Chi-Thanh Tran¹and Pavel Kudinov¹

Academic Editor: Xu Cheng

Received03 Apr 2013

Accepted07 May 2013

Published30 May 2013

Abstract

This paper is concerned with the development of approaches for assessment of core debris heat transfer and Control Rod Guide Tube (CRGT) cooling effectiveness in case of a Boiling Water Reactor (BWR) severe accident. We consider a hypothetical scenario with stratified (metal layer atop) melt pool in the lower plenum. Effective Convectivity Model (ECM) and Phase-Change ECM (PECM) are developed for the modeling of molten metal layer heat transfer. The PECM model takes into account reduced convection heat transfer in mushy zone and compositional convection that enables simulations of noneutectic binary mixture solidification and melting. The ECM and PECM are (i) validated against relevant experiments for both eutectic and noneutectic mixtures and (ii) benchmarked against CFD-generated data including the local heat transfer characteristics. The PECM is then applied to the analysis of heat transfer in a stratified heterogeneous debris pool taking into account CRGT cooling. The PECM simulation results show apparent efficacy of the CRGT cooling which can be utilized as Severe Accident Management (SAM) measure to protect the vessel wall from focusing effect caused by metallic layer.

1. Introduction

We consider a hypothetical severe accident in a BWR with subsequent core degradation, melt relocation, and debris bed (or cake) formation in the lower plenum filled with water. In case of inadequate cooling the debris bed (cake) will be heated up and remelted and a melt pool(s) will be formed. Prediction of transient melt pool formation, thermomechanical loading on the vessel and subsequent vessel failure modes [1], and timing and melt discharge characteristics is of paramount importance for the ex-vessel melt risk quantification in the Swedish BWR with a deep water-filled cavity under the reactor [2].

The lower plenum of a BWR contains a forest of CRGTs. In normal operation there is a purging water flow into the reactor through the CRGTs. In a severe accident the CRGT purging flow can be used for cooling the core melt materials and thus to become a potentially effective SAM measure for Swedish BWRs. Namely, the CRGT cooling may help to remove effectively the decay heat from a debris bed or melt pool formed in the lower plenum and thus delay or even prevent vessel failure [3] leading in the last case to in-vessel melt retention.

Besides the CRGT cooling, the other factors which may affect in-vessel progression of an accident are the phase changes involved in the melt pool formation process, and pool stratification (with separation of oxidic and metallic layers). There are large aleatory and epistemic uncertainties in the scenario and phenomena of melt pool formation, especially in the presence of core material physicochemical interactions [4]. However, homogeneous and horizontally stratified melt pools are the two most probable configurations which are considered in the present work.

In case of melt pool formation in the lower plenum, direct simulation of flow and heat transfer is a computational challenge, due to large length scale of the reactor lower plenum and high Rayleigh numbers (10¹⁵–10¹⁷ of the oxidic pool and 10⁸–10¹¹ of the metal layer). The difficulty is augmented with the presence of phase changes, complex 3D geometry of the lower head with a forest of CRGTs, complex flow patterns, and heat transfer induced by the CRGT cooling.

For prediction of decay-heated homogeneous melt pool formation in the Light Water Reactor (LWR) lower plenum, recently the ECM and PECM have been developed [5–7]. However, the stratified melt pool configuration with a metal layer atop which is characterized by Rayleigh-Benard convection or mixed natural convection (i.e., in the presence of side cooling, e.g., CRGT cooling) has not been considered. It is necessary to extend the ECM method to a metal layer. In the present work we develop prediction tools which enable simulations of the melt pool formation and melt-structure-water interaction with taking into account the CRGT cooling, homogeneity and stratification of melt pool configurations, and phase changes during the melt pool formation process.

The technical approach adopted in the present study is as follows. ECM and PECM are further developed to enable simulations of metal layer heat transfer on the base of effective convectivity approach [5, 6]. The CFD study is performed to gain insights into flow physics and examine Rayleigh-Benard convection and mixed natural convection heat transfer. The CFD study on the one hand, supports selection of appropriate correlations to be implemented in the ECM/PECM. On the other hand, CFD method is used to generate spatial distribution data of local heat transfer characteristics which are necessary for validation of the ECM and are not recoverable by average heat transfer correlations. The metal layer ECM and PECM are then validated against heat transfer experiments for both eutectic and non-eutectic binary mixture (melt), as well as against the CFD-generated data. The validated PECM is used to simulate heat transfer of a reactor-scale stratified melt pool whose formation is probable upon a certain BWR severe accident scenario which is also discussed in the paper.

The structure of the paper is as follows. In Section 2, we introduce development of the metal layer ECM and PECM. In Section 3 the CFD study is presented. Validation of the metal layer ECM and PECM is shown in Section 4. Section 5 contains the PECM application and discussions. Section 6 provides a brief summary of the paper.

2. Development of the Metal Layer Effective Convectivity Model

2.1. Review of Available Methods for Prediction of Metal Layer Heat Transfer

Metal layer heat transfer can be predicted using several methods which can be grouped into two classes: the lumped parameter and distributed parameter methods.

The lumped-parameter method is widely used [8–10] for calculation of heat fluxes from a metal layer to surroundings. In the lumped parameter method, the heat fluxes at the metal layer boundaries are calculated using the energy balance equation and empirical heat transfer correlations.

To determine the upward heat transfer coefficient (through the layer), the Globe-Dropkin correlation [11] is used in the work of Theofanous et al. [8], SCDAP/RELAP5-3D [9]. In MELCOR [10] code the Globe-Dropkin correlation is applied for the lower surface and the correlation produced in the ACOPO experiments [12] is applied for the upper surface of the metal layer. Modified Globe-Dropkin correlation which accounts for conductive heat transfer by adding the unity to the Nusselt number is employed in the MAAP code [13].

For determination of the sideward heat transfer coefficients a wider spectrum of different correlations is used. In the work of Theofanous et al. [8] and SCDAP/RELAP5-3D [9], the sideward heat transfer coefficient along a vertical cooled wall is the average Churchill-Chu correlation [14], while in the MELCOR code [10] it is the ACOPO correlation [12]. In the MAAP code, the sideward heat transfer coefficient is a function of the heat flux transferred at the lower boundary of the oxidic pool to the vessel wall [13]. In a later study which focused on the implemented correlations in MAAP applied to French PWRs [15], the sideward heat transfer model was replaced by the Churchill-Chu correlation which results in higher probability of focusing effect in MAAP calculations for French PWRs.

Clearly, the lumped-parameter method is a simple and effective method for calculations of heat fluxes and energy splitting. Note that the lumped-parameter method established in [8] was validated against the MELAD experiments which were built specially for validation of mixed natural convection heat transfer models. However, the phase change crust dynamics are not considered in the lumped-parameter methods. Thus the effect of latent heat of fusion on the transient heat fluxes at the metal layer boundaries cannot be quantified (we will show the effect of latent heat on heat fluxes in this paper). Furthermore, the local effect of the heat flux along an inclined cooled surface is not considered in the method.

The CFD method such as Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) is extensively used to examine flow characteristics; low-Prandtl number effects of a molten metal fluid layer heated from below and cooled from the top [16–19]. To understand complex fluid flows, the CFD method is indispensable and helps to reveal separate-effect phenomena. However, the CFD method is computationally expensive and thus difficult to be applied to simulations of melt pool heat transfer. The difficulty is further increased for long (hours) transient scenarios, for the cases of parametric studies, and especially when the phase change is involved.

The Effective Convectivity Conductivity Model (ECCM) can be classified as a distributed parameter method and is a computationally affordable alternative of CFD simulations. Originally the ECCM was developed by Bui and Dinh [20], later on extended to simulations of metal layer heat transfer, and reported in detail in the work of Sehgal et al. [21]. The metal layer ECCM employs the effective convectivity model to represent turbulent heat transport (Rayleigh-Benard natural convection) through a fluid layer using so-called effective velocity which is defined as follows: where is the thickness of the fluid layer and is empirical Rayleigh-Benard convection heat transfer coefficient. To describe turbulent heat transport in the horizontal direction, the effective conductivity model is used (). In contrast to the lumped-parameter methods, the boundary layer model [22] is applied to describe the local effect of the heat transfer coefficient along a vertical/inclined cooled wall:

The metal layer ECCM was implemented in MVITA code which was applied to simulations of stratified melt pool heat transfer in a PWR lower head and to examine the focusing effect. However, MVITA is a 2D code, thus it is not applicable to 3D geometry of a BWR lower plenum. Furthermore, in the ECCM non-eutectic mixture phase change and mushy zone heat transfer were not considered.

2.2. Governing Equations and Heat Transfer Characteristic Velocities

Following the concept of effective convectivity [6, 20] the ECM is developed to enable effective heat transfer simulations for a metal layer which may appear atop of an oxidic melt pool in the lower plenum.

The metal layer ECM uses the directional heat transfer characteristic velocities to describe turbulent heat transfer in both the upward and sideward (in the presence of side cooling) directions of a fluid layer. As the heat transfer characteristic velocities (later on denoted as and ) are used instead of instantaneous velocity hence, the following energy conservation equation is solved:

The heat transfer characteristic velocities are determined as follows (see Figure 1).

For a fluid layer heated from below and cooled from the top, using the upward heat transfer characteristic velocity , the amount of convective heat transferred from the bottom surface to the bulk fluid (and from the bulk fluid to the top surface) is defined as follows: where is the area of the cooled top/heated bottom surface and is the driving temperature difference between the bottom and top surfaces.

Conduction heat transfer from the bulk fluid to the top cooled surface is

The heat balance equation is then written as follows:

From (6), the upward heat transfer characteristic velocity is derived as follows:

The heat transfer coefficient in (7) is defined using external Rayleigh number Ra which is a function of the temperature difference across the layer.

In the case of mixed cooling by the top wall and one side wall (other side wall is either adiabatic or symmetrical, Figure 1), using a formal analogy, the sideward characteristic velocity is derived as follows:

Note that the heat transfer coefficient is defined using the sideward Rayleigh number which is a function of the sideward temperature drop, that is, the difference between the bulk temperature and vertical/inclined cooled boundary temperature. Selection of the upward and sideward heat transfer correlations and description of the sideward heat transfer coefficient profile are based on the findings and insights gained from the CFD study which is presented in Section 3.

To implement the ECM, the convective terms with the heat transfer characteristic velocities in (3) are defined as a modified heat source term. In such a way (3) is reduced to a heat conduction equation, can be solved using a commercial CFD code, for example, the Fluent code.

For a fluid layer heated from below and cooled from the top (and side) walls, the characteristic velocity is positive on a cooled surface and negative on a heated surface (Figure 1). Apparently the convective terms added to the grid cells, where temperature gradients are not zero, are artificial. In order to keep energy balance of the computational domain, the fluid bulk temperature must be adjusted based on the summarized heat of the artificial convective terms which have been added to the grid cells in the previous time step.

2.3. Phase-Change Effective Convectivity Model

In the phase-change ECM (PECM) for a metal layer, the single enthalpy conservation equation which is common for solid, mushy, and liquid phases is solved:

Similarly to the PECM for an internally heated volume [7], (9) is solved using the source-based method with a fixed grid [23]. The main approach for the metal layer treatment in PECM is to define the heat transfer characteristic velocities in a mushy zone () and to employ the compositional convection model.

Heat transfer characteristics of a mushy zone depend on the certain fluid velocity in it. One may assume that the fluid velocity is a function of liquid fraction which in turn is related to temperature of the binary mixture in a mushy zone by either linear, or linear-eutectic, or Scheil, or power-law relationship [23]. The liquid fraction can be determined as follows:

Thus in the PECM, based on (10) different models of the mushy zone characteristic velocities can be realized using different types of liquid fraction dependency, that is, . Therefore in a mushy zone, reduced characteristic velocities are employed to describe mushy zone heat transfer.

Natural convection in a fluid layer involving non-eutectic mixture phase change is affected not only by the temperature difference but also by the compositional concentration difference of a solute across the liquid layer (double-diffusive convection). In such a case, the effective Rayleigh number which defines the intensity of turbulent natural convection in a fluid layer can be determined as a combination of the thermal Rayleigh number Ra_T and compositional Rayleigh number Ra_C. Thermal Rayleigh number Ra_T is defined as follows [24]: where is defined as . Compositional Rayleigh number Ra_C can be determined as: where is the concentration difference across the fluid (liquid) layer and is the rejectability coefficient. The presence of the rejectability coefficient in the compositional Rayleigh number expression is explained as follows.

According to the definition, the concentration difference is . However, binary mixture solidification experiments [25, 26] show that the value of remains uncertain due to rejectability of the higher concentration (compared with the bulk fluid) solute from a mushy zone to the bulk fluid. For one binary mixture, the solute with higher concentration may easily be rejected from the mushy zone, while for the other mixture, this higher concentration solute may be stuck in the mushy zone and be solidified. The rejectability of the solute from a mushy zone also depends on the mushy zone structure, the intensity of natural convection in both the mushy zone and bulk fluid, and the fluid properties and the cooling rate. Note that depending on the hypereutectic or hypoeutectic mixture and the cooling direction, that is, cooling from the upper or lower surface, compositional convection may weaken/strengthen thermal convection. We envision the necessity of examining the effect of mushy zone heat transfer characteristic velocity models and rejectability coefficient on heat transfer of non-eutectic binary mixture solidification/melting.

In the present paper, we show that with appropriate selection of the liquid fraction dependence model for the mushy characteristic velocity, and the rejectability coefficient, the PECM is able to predict the transient behavior of non-eutectic binary mixture solidification observed in different experiments (Section 4).

3. CFD Study

The main purpose of the CFD study is to examine Rayleigh-Benard convection and mixed natural convection heat transfer in order to support selection of appropriate correlations to be implemented in the metal layer ECM. Also CFD data help to quantify the local sideward heat transfer coefficient distribution for validation of the ECM.

The CFD method used in the present study is named Implicit LES (ILES) [27]. The ILES method uses a high-resolution grid to effectively provide LES without explicit Subgrid Scale (SGS) turbulence modeling. To capture the wall-boundary effects, a higher-resolution grid is provided in the near-wall region of the computational domain. The results obtained with ILES for heat transfer coefficients are grid independent.

Prior to presenting CFD simulation results, let us make a brief summary of the past works related to Rayleigh-Benard turbulent convection with respect to low Pr number fluids. It was found that heat transfer coefficient depends on Rayleigh number (external) and the dependency is an exponential law function, where the value of the exponent ranges from 1/4 to 1/3 [28]. In the work of Goldstein et al. [29], it was predicted that the experimental exponent of Ra seems to increase with increasing Pr number. Very weak dependence of the Nusselt number on Pr is observed in the range of [30]. However, for low Pr number fluids Nusselt number becomes a function of not only Ra number but also of Pr number.

In the present study, the CFD predicted results are compared with three correlations which were obtained in the experiments for low Prandtl number fluids. The first is Globe-Dropkin correlation [11], the second is O’Toole-Silveston correlation [31], and the third is Cioni correlation [32]. The Globe-Dropkin correlation is given as follows:

The O’Toole-Silveston correlation is

The Cioni correlation is expressed as follows:

The CFD method (ILES) is applied to simulations of low Pr number fluid layer (liquid metal) in rectangular cavity and Unit Volume (UV) geometries. The UV is a representative rectangular cavity of molten metal surrounding a CRGT. The height of UV is varying depending on the metal layer thickness. To examine the effect of aspect ratio (the height-to-width ratio), the CFD simulations are performed for different heights of rectangular cavity and UV.

Figure 2 shows plots of heat transfer coefficients across the layer heated from below and cooled from the top (Rayleigh-Benard convection), calculated by heat transfer correlations (13), (14), and (15) and the CFD predicted data for different geometries. The predicted Nusselt number is well agreed with the Cioni correlation for the rectangular cavity with the aspect ratio of less than unity (Figure 2). For the UV geometry, two predicted Nusselt numbers shown in the figure are for two different aspect ratios. In the UV with an aspect ratio larger than unity (the upper circle dot), the Nusselt number is consistent with the Globe-Dropkin correlation, while in the other (the lower circle dot with an aspect ratio of less than unity) the Nusselt is slightly lower and consistent with the Cioni correlation. Apparently the aspect ratio does affect the heat transfer; however, the influence can be assumed insignificant.

To examine the effect of side cooling on heat transfer coefficients, a CFD simulation is performed for a fluid layer cooled by both the top and side walls (mixed natural convection). Figure 3 presents temperature profiles across the metal layer of two cooled configurations: the first is Rayleigh-Benard convection and the second is mixed natural convection (i.e., with side cooling) in the same geometry. Clearly the profiles are qualitatively similar (Figure 2). In case of side cooling, the bulk temperature of the layer is decreased that results in changing the heat transfer coefficients along the top and bottom walls.

Figure 4 presents the CFD predicted heat transfer coefficients of a mixed natural convection fluid layer. As it is shown in the figure, the heat transfer coefficients along the upper surface are lower than those of the lower surface (the heated surface). The CFD predicted Nusselt numbers are fairly consistent with the experimental correlations. It is seen that the heat transfer coefficients predicted for the UV are higher than those of the rectangular cavity and consistent with the Globe-Dropkin correlations. This may be related to specific geometry of the UV, that is, the presence of a CRGT which changes flow structure and thus results in more intensive heat transfer.

For the mixed natural convection fluid layer, the sideward heat transfer coefficients (along the vertical wall and UV cooled CRGT wall) are presented in Figure 5. The predicted Nusselt numbers are well agreed with the Churchill-Chu correlation.

Figure 6 shows the heat flux profile along the CRGT (mixed natural convection in UV geometry) and the analytical model which is derived from the average Churchill-Chu correlation [14] and the boundary layer model (2). The two profiles are found to be in good agreement along the vertical wall. Although in the uppermost and lowermost regions, the heat flux values are slightly divergent. An explanation of these deviations is as follows. At the upper part of CRGT surface the heat flux is decreased significantly due to the isothermal boundary conditions (with the same temperature) applied for both top (horizontal) and CRGT (vertical) walls. For the lower part of the vertical wall, due to the cold fluid descending from the boundary layer, temperature of the mixed fluid flow is decreased, causing decrease of the heat flux along this part (smaller temperature difference). Note that rapid increase of the heat flux obtained with the analytical model at the top of CRGT is artificial. However, the analytical model presents a conservative estimation of local heat flux.

On the basis of analysis of the CFD simulations results we adopt the Globe-Dropkin correlation for the upward heat transfer coefficient and the Churchill-Chu correlation for the sideward heat transfer coefficient in the ECM and PECM. We select these correlations because they are consistent with the CFD simulation results obtained for the specific BWR geometries under mixed natural convection conditions with low Pr number fluids. Furthermore, these correlations are valid in a wide range of Rayleigh number which is close to the prototypic reactor condition and applicable for low Pr number fluids (e.g., liquid metals). The profile of the sideward heat transfer coefficient can be described by the boundary layer model, (2) which is also confirmed by the CFD study with geometry of interest.

For validation of the developed method against phase change heat transfer we consider simulation of an experiment on solidification and heat transfer in a fluid layer [33]. The experiment was performed in a rectangular cavity of 6.35 cm × 3.81 cm × 8.89 cm; the simulant used was pure gallium (Ga). The cavity is cooled from the top surface till freezing while the bottom surface is heated at constant temperature. Temperature profiles across the layer and crust thickness evolution during the experiment were reported.

The CFD ILES method is applied to simulate this experiment. Figure 7 shows the CFD predicted velocity and liquid fraction contours in 22 min. A Benard’s cell is clearly shown in the figure. Temperature profiles across the layer (in 2 min and 8 min) are presented in Figure 8, along with the experimental data. Apparently, the CFD predicted results are well agreed with the experimental data. The CFD simulation data are used for validation of the PECM (see the next section).

(a)

(b)

4. Validation of the ECM and PECM

Validation of the ECM and PECM covers a wide spectrum of physical phenomena involved in metal layer heat transfer such as Rayleigh-Benard convection mixed natural convection boundary layer development and transient phase change and crust formation. As the dual approach is adopted in the present study, both experimental and CFD-generated (with DNS and ILES methods) data are used for the validation of the ECM and PECM models. Various heat transfer problems are considered in the validation matrix (Table 1), including classical Rayleigh-Benard convection, mixed natural convection and as well as eutectic and non-eutectic binary mixture solidification.

4.1. ECM Validation

DNS was used in [34] to predict the temperature profile of a fluid layer heated from below and cooled from the top. The fluid Pr number was 0.022 and the layer Rayleigh number reached 2.2 × 10⁷. The metal layer ECM is used to predict the temperature profile across the fluid layer. Figure 9 presents temperature profiles predicted by the ECM and DNS methods. It can be seen in the figure, good agreement of the two predicted results is obtained.

The ECM is also used to predict the Rayleigh-Benard convection temperature profile with a lower Ra number (Ra = 6.3 × 10⁵) and of order of unity Pr number (). Figure 10 shows that the ECM predicted the temperature profile agrees well with that predicted by the DNS method [19].

To validate the ECM in prediction of temperature and energy splitting of a mixed natural convection fluid layer, the MELAD A1 experiment [8] is used. The ECM predicted results are presented in Table 2. Clearly, the ECM predicted temperature and heat fluxes are very close to the experimental ones. A slight difference in the bulk temperature is observed. This may relate to possible slight thermal stratification in the fluid layer which is not captured by the ECM method.

As the next step, the ECM is benchmarked against CFD simulation data. It is well known that CFD methods can be used to perform reliable “numerical experiments” for classical Rayleigh-Benard convection [16]. In the present study, the CFD method (ILES) is used to simulate heat transfer of a mixed natural convection fluid layer. A UV of 31.26 cm × 31.26 cm × 40 cm is heated from below and cooled from the top wall and from the CRGT wall. The obtained steady state temperature profile and heat flux distribution along the CRGT wall are used for ECM validation.

Figure 11 presents plots of the CFD and ECM predicted temperature profiles across the layer. The two profiles are very close across the most fluid layer. Figure 12 shows the heat flux profiles along the CRGT wall; they are well agreed in the middle part of the CRGT wall. In the upper region, the ECM predicted heat flux value is slightly higher than that obtained with the CFD (less than 10%). In the lower region of the cooled wall, the CFD method is able to predict diminishing heat flux due to the flow stagnation, while the ECM is unable to capture this effect. However, ECM estimation of heat flux can be considered as reasonable and conservative in the sense that the heat flux is not underestimated by ECM model.

4.2. PECM Validation

In this section, validation of the PECM against both eutectic and non-eutectic binary melts is presented. First, the PECM is used to simulate a eutectic melt solidification experiment [33] which was also simulated by the CFD ILES method.

Figure 13 presents crust thickness evolutions predicted by the PECM, CFD, and the experimental data. Excellent agreement between the methods is obtained. In the later time period (after 20 min), the CFD predicted crust is slightly thinner than that of the experiment and the PECM simulation. This is the result of a larger CFD predicted heat transfer coefficient across the convective fluid layer. However, the difference between the CFD and experimental data is insignificant.

The PECM is also benchmarked against the CFD-generated data (Section 3). Figure 14 shows the PECM simulation temperature profile along with that of the CFD simulation. The temperature profiles predicted by the two methods are in good agreement. Note that the CFD simulation of the experiment takes about 5 days of 2.8 GHz CPU time while the PECM simulation lasts only few hours.

To validate the mushy characteristic velocity model and rejectability coefficient implemented in the PECM we use data from one of the experiments of Cao and Poulikakos [35]. In the experiments, non-eutectic binary mixtures of the aqueous ammonium chloride (NH₄Cl–H₂O) in a cavity of 48.3 cm × 25.4 cm × 12.7 cm are cooled from the top surface till freezing. Different compositions of ammonium chloride were used. Evolutions of crust thickness (solid and mushy) were recorded.

The PECM is used to simulate the experiment where the NH₄Cl composition is 5%. PECM simulations show that crust thickness evolution is sensitive to different models of mushy characteristic velocity and to different values of rejectability coefficient. Figure 15 presents the experimental crust thickness evolution and results predicted by the PECM simulation where the mushy characteristic velocity model is described by a third-order power law of liquid fraction and the rejectability coefficient is 0.07. Although there is a slight deviation in prediction of the solid thickness start point (in around 4 h), the PECM well captures the evolutions of the both solid and mushy crusts thicknesses.

Extensive validation of the ECM and PECM against various experiments and CFD-generated data has confirmed the applicability of the ECM and PECM for simulations of turbulent natural convection heat transfer in a fluid layer heated from below, cooled from the top and side walls, as well as for simulations of transient heat transfer which involves the phase change in a BWR lower plenum configuration.

5. Application of the PECM to BWR Melt Pool Simulations

5.1. A Severe Accident Scenario in a BWR

An accident scenario which can lead to formation of stratified melt pool (Figure 16) in the BWR lower head is considered in this section. The scenario is started from state 1, when a debris bed (cake) is formed in the water-filled lower plenum. It is assumed that water ingress is possible only in the upper layer of the debris bed, resulting in formation of a coolable debris layer (bed) atop, while the lower part of the bed is reheated up with time. As temperature of the bed (lower part) exceeds liquidus temperature of metals, metallic components of the bed will be melted and flowing down, filling the pores in the lower part of the debris bed. Due to much higher liquidus temperature the ceramic component will remain in the solid phase (particulates). It can be assumed that displacement of solid oxidic particulates will be possible due to the gradual motion downwards and avalanches. Consequently a heterogeneous debris pool (i.e., oxidic particulates in a molten metal pool) will be formed in the lower plenum (state 2 in Figure 16). Depending on the mass fraction of metallic components in the quenched debris, the debris pool may become stratified; that is, oxidic particulates will be submerged in a liquid metal pool with a metal layer floating atop. In case of insufficient molten metals to fill the pores in the debris bed, a uniform (unstratified) heterogeneous debris pool will be formed in the lower plenum.

In the present study, we apply the metal layer PECM to simulate heat transfer of the stratified heterogeneous debris pool. It is assumed that the IGTs are intact or plugged and do not affect heat transfer of the pool, and chemical interaction is not considered.

5.2. PECM Simulation, Results, and Discussion

The PECM heat transfer simulation of stratified heterogeneous debris pool is presented in this section. The stratified heterogeneous debris pool is comprised of a debris pool which is characterized by a matrix of solid particulates and molten metals and a molten metal layer atop of the pool. The metal layer thickness depends on the fraction of the metallic components presented in the initial debris bed with the porosity of 40% (Figure 17). As an example, we consider the debris pool depth of 0.8 m and the metal layer thickness of 0.2 m (total m). The decay heat (about 1 MW/m³) remains only in the debris pool. The simulation is performed for a representative slice of the BWR lower plenum. The slice is a 3D segment which is filled with core materials (i.e., solid ceramic particulates and liquid metals) and includes 6 cooled CRGTs and a section of the vessel wall from below (Figure 18).

The simulation is started from a hot condition when temperature of the debris pool is higher than liquidus temperature of metals (2000 K is applied). In the presence of CRGT cooling, it is supposed that initial crusts are available on the CRGT wall surfaces. The thickness of the initial crusts is predicted to be of 40 mm, based on the simulation performed for a uniform heterogeneous debris pool. The remained volume of the metal layer is hot and in liquid state. Such a scenario is rather of low probability because formation of metal layer is a process which takes time. However, we consider such a “hot state” as a bounding scenario of instantaneous metal layer formation. Due to low permeability of the porous media, it is supposed that there is no convection in the debris pool; only heat conduction is considered. Conductivity of the debris pool is defined as the effective conductivity of the ceramic oxide and metals (). Natural convection is assumed to take place only in the molten volume of the metal layer where the PECM is applied.

The boundary conditions are as follows. For the CRGT wall internal surfaces the isothermal boundary condition is applied. The top surface of the metal layer is in close contact with the decay heated coolable debris layer (Figure 16), the lower boundary of which may be heated up to high temperature (e.g., liquidus temperature of metals). Therefore the radiation heat transfer boundary condition is used for the metal layer top surface. External radiation temperature is assumed to be liquidus temperature of metals for the sake of conservatism. The external surface of the vessel wall is assumed to be insulated during the accident, thus a small heat flux (about 50 W/m²) is allowed on this surface. The other boundaries are adiabatic or symmetrical.

Figure 18 presents a transient state of the upper metal layer and the debris pool. There are growing crusts on CRGT surfaces and on the vessel wall submerged in the metal layer. Note that the debris pool state remains unchanged, that is, a matrix of particulates and molten metals, although it is shown in the figure as solid.

The PECM simulation shows that after about 30 min, although temperature of the lower debris pool (core) becomes higher than 2100 K, the upper metal layer is fully solidified. This is because of efficient heat removal from the high thermal conductivity metal layer to the cooled CRGTs. Steady state temperature distribution in the pool is presented in Figure 19 where the metal layer is shown at low temperature. Steady state temperature of the lower debris pool is about 2040 K indicating that the metal components remain in the liquid phase, and ceramic oxidic particulates are not melted.

Figure 20 presents steady state temperature profiles of the vessel wall external surfaces for two cases under the same total amount of heat generation and pool depth (1.0 m): the first is the considered stratified heterogeneous debris pool ( MW/m³ in the lower part and MW/m³ in the metal layer), and the second is the uniform heterogeneous debris pool ( MW/m³), without a metal layer atop. In the first case, due to effective CRGT cooling, the heat flux to the vessel is low and vessel wall temperature is kept well below 1100°C, while in the second case, even with a smaller volumetric heat generation rate, temperature of the vessel approaches 1200°C. Creep temperature for the reactor vessel steel (e.g., SA533B1) usually assumed at 1100°C [1]. Thus vessel failure is not predicted upon the stratified heterogeneous debris pool scenario. However, for the uniform heterogeneous debris pool, the vessel wall may fail in the place connected to the uppermost region of the debris pool.

Note that the CRGT cooling is provided by the inside water flow; thus it is efficient only when the heat flux to CRGTs is lower than the Critical Heat Flux (CHF) which in turn depends also on water flow rate. It is therefore important to identify whether the heat flux along a CRGT is higher than the CHF.

The PECM simulation shows that, for the stratified heterogeneous debris pool scenario, the steady state heat flux to the CRGTs (234 kW/m²) is lower than the CHF at nominal CRGT purging flow rate of 15 kg/(m²·sec) (around 400 kW/m²). Due to the CRGT cooling efficacy and due to the low heat flux to the metal layer from the below debris pool, the focusing effect is not observed.

Heat flux evolution (Figure 21) shows that the transient heat flux to the CRGT is highest during the first 20 min. The source of the high heat flux is the latent heat released during fast solidification of metals (Figure 21). However, this transient heat flux to the CRGTs is well below the CHF for the CRGTs at water flow rate of 30 kg/(m²·sec) (about 900 kW/m²) [36]. This means the effect of latent heat is small as long as we can provide a sufficient water flow rate to avoid CHF (about ~2–4-folds of nominal water flow rate).

It is worth to mention that at the nominal purging flow rate of ~15 kg/(m²·sec) when the CHF can be considerably smaller [36], failure of the cooled CRGT upon the transient heat flux remains questionable. For more rigorous analysis we envision a need to examine the effect of variable boundary conditions on the transient heat flux to CRGT walls. The study approach is coupling simulation [37] using the metal layer PECM and RELAP which is applied for CRGTs with varying water follow rates to 2–4-folds of the nominal one.

6. Concluding Remarks

The metal layer Effective Convectivity Model (ECM) and Phase-change ECM (PECM) presented in this paper are built on a concept of the previously developed Effective Convectivity Conductivity Model (ECCM) [20] and follows the technical approach of recently developed ECM and PECM for heat transfer simulations of a decay-heated melt pool [6, 7].

Insights gained from the Computational Fluid Dynamics (CFD) study helped to improve the metal layer ECM and PECM. The metal layer PECM is augmented by the mushy zone heat transfer characteristic velocity and compositional convection models for simulation of non-eutectic binary mixtures.

Validation of the developed models is performed against both experimental and CFD-generated data. It is demonstrated that the modified metal layer ECM and PECM are capable of predicting energy splitting in metal layers for various heat transfer problems with reasonable accuracy.

The PECM application to heat transfer simulation in a BWR severe accident scenario shows the high potential of CRGT flow cooling as an efficient SAM measure to remove the decay heat from a stratified heterogeneous debris pool, keeping lower head vessel wall temperature well below the thermal creep limit. However, failure of the vessel should be considered using coupled thermo-mechanical creep analysis [38–40]. Further analysis will be performed for assessment of necessary CRGT flow rate which will ensure with sufficient margin no failure of CRGTs.

Nomenclature

Arabic

	Composition, %
:	Specific heat capacity, J/(kg·K)
:	Liquid fraction
:	Gravitational acceleration, m/s²
:	Depth (or thickness) of a fluid layer, m
:	Sensible enthalpy, J/kg
:	Conductivity, W/(m·K)
:	Nusselt number,
Pr:	Prandtl number,
:	Heat, J
:	Volumetric heat generation rate, W/m³
Ra:	External Rayleigh number,
:	Local Rayleigh number,
:	Rejectability coefficient
:	Area, m²
:	Temperature, K
:	Fluid velocity, m/s
:	Characteristic velocity, m/s
:	Mushy characteristic velocity, m/s
:	Width of a fluid layer, m
:	Local vertical coordinate, m.

Greek

Δ:	Latent heat, J/kg
Δ:	Temperature difference, K
:	Thermal diffusivity, m²/s,
:	Thermal expansion coefficient, 1/K
:	Kinematics viscosity, m²/s
:	Density, kg/m³
:	Bulk.

Subscripts and Superscripts

:	Compositional
cond:	Conduction
conv:	Convection
:	Eutectic
eff:	Effective
:	Componentindex
LIQ:	Liquidus
RB:	Rayleigh-Benard
side:	Sideward
SOL:	Solidus
:	Thermal
up:	Upward
:	Axis direction.

References

J. L. Rempe, S. A. Chavez, G. L. Thinnes et al., “Light water reactor lower head failure analysis,” Tech. Rep. NUREG/CR-5642, EGG-2618, Idaho National Engineering Laboratory, Washington, DC, USA, 1993.
View at: Google Scholar
P. Kudinov, A. Karbojian, W. M. Ma, and T. N. Dinh, “An experimental study on debris formation with corium stimulant materials,” in Proceedings of the International Congress on Advances in Nuclear Power Plants (ICAPP '08), Anaheim, Calif, USA, June 2008.
View at: Google Scholar
C. T. Tran and T. N. Dinh, “Application of the phase-change effective convectivity model to analysis of core melt pool formation and heat transfer in a BWR lower head,” in Proceedings of the Annual Meeting of the American Nuclear Society, pp. 617–618, Anaheim, Calif, USA, June 2008.
View at: Google Scholar
V. G. Asmolov, S. V. Bechta, V. B. Khabensky et al., “Partitioning of U, Zr and Fe between molten oxidic and metallic corium,” in Proceedings of the MASCA Seminar, Aix-en-Provence, France, June 2004.
View at: Google Scholar
C. T. Tran, P. Kudinov, and T. N. Dinh, “An approach to numerical simulation and analysis of molten corium coolability in a boiling water reactor lower head,” Nuclear Engineering and Design, vol. 240, no. 9, pp. 2148–2159, 2010.
View at: Publisher Site | Google Scholar
C. T. Tran and T. N. Dinh, “An effective convectivity model for simulation of in-vessel core melt progression in boiling water reactor,” in Proceedings of the International Congress on Advances in Nuclear Power Plants (ICAPP '07), Nice Acropolis, France, May 2007.
View at: Google Scholar
C. T. Tran and T. N. Dinh, “Simulation of core melt pool formation in a Reactor pressure vessel lower head using an Effective Convectivity Model,” Nuclear Engineering and Technology, vol. 41, no. 7, pp. 929–944, 2009.
View at: Google Scholar
T. G. Theofanous, C. Liu, S. Additon, S. Angelini, O. Kymalainen, and T. Salmassi, “In-vessel coolability and retention of a core melt,” DOE/ID-1046, 1994.
View at: Google Scholar
SCDAP/RELAP5-3D Code Development Team, “SCDAP/RELAP5-3D Code Manual,” Report INEEL/EXT-02-00589, Revision 2.2, Idaho National Engineering and Environmental Laboratory, 2003.
View at: Google Scholar
R. O. Gauntt, R. K. Cole, C. M. Erickson et al., MELCOR Computer Code Manual, Core (COR) Package Reference Manuals, NUREG/CR-6119, vol. 2, Rev. 2, Version 1.8.6, 2005.
S. Globe and D. Dropkin, “Natural-convection heat transfer in liquid confined by two horizontal plates and heated from below,” Journal of Heat Transfer, vol. 81, pp. 24–28, 1959.
View at: Google Scholar
T. G. Theofanous, M. Maguire, S. Angelini, and T. Salmassi, “The first results from the ACOPO experiment,” Nuclear Engineering and Design, vol. 169, no. 1–3, pp. 49–57, 1997.
View at: Google Scholar
MAAP4 Users Manual, vol. 2, Fauske Associated, 1999.
S. W. Churchill and H. H. S. Chu, “Correlating equations for laminar and turbulent free convection from a vertical plate,” International Journal of Heat and Mass Transfer, vol. 18, no. 11, pp. 1323–1329, 1975.
View at: Google Scholar
M. Eddi, “Study on heat transfer in lower head of nuclear power plant vessel during a severe accident,” in Proceedings of 9th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH '99), San Francisco, Calif, USA, October 1999.
View at: Google Scholar
G. Grotzbach, “Direct numerical simulation of laminar and turbulent benard convection,” Journal of Fluid Mechanics, vol. 119, pp. 27–53, 1982.
View at: Google Scholar
R. Verzicco and R. Camussi, “Prandtl number effects in convective turbulence,” Journal of Fluid Mechanics, vol. 383, pp. 55–73, 1999.
View at: Google Scholar
R. M. Kerr and J. R. Herring, “Prandtl number dependence of Nusselt number in direct numerical simulations,” Journal of Fluid Mechanics, vol. 419, pp. 325–344, 2000.
View at: Google Scholar
I. Otic, G. Grotzbach, and M. Worner, “Analysis and modelling of the temperature variance equation in turbulent natural convection for low-prandtl fluids,” Journal of Fluid Mechanics, vol. 525, pp. 237–261, 2005.
View at: Google Scholar
V. A. Bui and T. N. Dinh, “Modeling of heat transfer in heated-generating liquid pools by an effective diffusivity-convectivity approach,” in Proceedings of 2nd European Thermal-Sciences Conference, pp. 1365–1372, Rome, Italy, 1996.
View at: Google Scholar
B. R. Sehgal, V. A. Bui, T. N. Dinh, and R. R. Nourgaliev, “Heat transfer process in reactor vessel lower plenum during a late phase of in-vessel core melt progression,” Advances in Nuclear Science and Technology, vol. 26, pp. 103–135, 1998.
View at: Google Scholar
E. R. G. Eckert and T. W. Jackson, “Analysis of turbulent free convection boundary layer on flat plate,” NACA Technical Note 2207, 1950.
View at: Google Scholar
V. R. Voller and C. R. Swaminathan, “General source-based method for solidification phase change,” Numerical Heat Transfer B, vol. 19, no. 2, pp. 175–189, 1991.
View at: Google Scholar
M. G. Worster, “Natural convection in a mushy layer,” Journal of Fluid Mechanics, vol. 224, pp. 335–359, 1991.
View at: Google Scholar
J. S. Wettlaufer, M. G. Worster, and H. E. Huppert, “Natural convection during solidification of an alloy from above with application to the evolution of sea ice,” Journal of Fluid Mechanics, vol. 344, pp. 291–316, 1997.
View at: Google Scholar
R. Trivedi, H. Miyahara, P. Mazumder, E. Simsek, and S. N. Tewari, “Directional solidification microstructures in diffusive and convective regimes,” Journal of Crystal Growth, vol. 222, no. 1-2, pp. 365–379, 2001.
View at: Publisher Site | Google Scholar
L. G. Margolin, W. J. Rider, and F. F. Grinstein, “Modeling turbulent flow with implicit LES,” Journal of Turbulence, vol. 7, pp. 1–27, 2006.
View at: Publisher Site | Google Scholar
S. Grossmann and D. Lohse, “Scaling in thermal convection: a unifying theory,” Journal of Fluid Mechanics, vol. 407, pp. 27–56, 2000.
View at: Google Scholar
R. J. Goldstein, H. D. Chiang, and D. L. See, “High-Rayleigh-number convection in a horizontal enclosure,” Journal of Fluid Mechanics, vol. 213, pp. 111–126, 1990.
View at: Google Scholar
P. E. Roche, B. Castaing, B. Chabaud, and B. Hebral, “Prandtl and rayleigh numbers dependences in rayleigh-benard convection,” Europhysics Letters, vol. 58, no. 5, pp. 693–698, 2002.
View at: Google Scholar
J. L. O'Toole and P. L. Silveston, “Correlations of convective heat transfer in confined horizontal layers,” AIChE Chemical Engineering Progress Symposium Series, vol. 57, no. 32, pp. 81–86, 1961.
View at: Google Scholar
S. Cioni, S. Ciliberto, and J. Sommeria, “Strongly turbulent Rayleigh-Bénard convection in mercury: comparison with results at moderate Prandtl number,” Journal of Fluid Mechanics, vol. 335, pp. 111–140, 1997.
View at: Google Scholar
C. Gau and R. Viskanta, “Effect of natural convection on solidification from above and melting from below of a pure metal,” International Journal of Heat and Mass Transfer, vol. 28, no. 3, pp. 573–587, 1985.
View at: Google Scholar
A. A. Mohamad and R. Viskanta, “Modeling of turbulent buoyant flow and heat transfer in liquid metals,” International Journal of Heat and Mass Transfer, vol. 36, no. 11, pp. 2815–2826, 1993.
View at: Google Scholar
W. Z. Cao and D. Poulikakos, “Solidification of an alloy in a cavity cooled through its top surface,” International Journal of Heat and Mass Transfer, vol. 33, no. 3, pp. 427–434, 1990.
View at: Google Scholar
C. T. Tran and T. N. Dinh, “Analysis of melt pool heat transfer in a BWR lower head,” in Transactions of ANS Winter Meeting, pp. 629–631, Albuquerque, NM, USA, November 2006.
View at: Google Scholar
F. Cadinu, C. T. Tran, and P. Kudinov, “Analysis of in-vessel coolability and retention with control rod guide tube cooling in boiling water reactors,” in Proceedings of the NEA/SARNET In-Vessel Coolability (IVC) Workshop, Issy-les-Moulineaux, France, October 2009.
View at: Google Scholar
W. Villanueva, C.-T. Tran, and P. Kudinov, “Coupled thermo-mechanical creep analysis for boiling water reactor pressure vessel lower head,” Nuclear Engineering and Design, vol. 249, pp. 146–153, 2012.
View at: Publisher Site | Google Scholar
C. Torregrosa, W. Villanueva, C.-T. Tran, and P. Kudinov, “Coupled 3D thermo-mechanical analysis of a nordic BWR vessel failure and timing,” in Proceedings of the 15th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH '13), Pisa, Italy, May 2013, Paper 495.
View at: Google Scholar
A. Goronovski, W. Villanueva, C.-T. Tran, and P. Kudinov, “The Effect of internal pressure and debris bed thermal properties on BWR vessel lower head failure and timing,” in Proceedings of the 15th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH '13), Pisa, Italy, May 2013, Paper 500.
View at: Google Scholar

Copyright

Copyright © 2013 Chi-Thanh Tran and Pavel Kudinov. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

2853

Downloads

1561

Citations