Abstract

This paper presents, the numerical analysis of heat and mass transfer during hydrogen generation in an array of fuel cylinder bars, each coated with a cladding and a steam current flowing outside the cylinders. The analysis considers the fuel element without mitigation effects. The system consists of a representative periodic unit cell where the initial and boundary-value problems for heat and mass transfer were solved. In this unit cell, we considered that a fuel element is coated by a cladding with steam surrounding it as a coolant. The numerical simulations allow describing the evolution of the temperature and concentration profiles inside the nuclear reactor and could be used as a basis for hybrid upscaling simulations.

1. Introduction

During a severe nuclear reactor accident, several works are focusing on core degradation by metal core components oxidation by air or steam [14]. Studies on possible consequences of core meltdown have demonstrated that hydrogen combustion is one of the contributors to the containment early failure [5]. Justly, during Fukushima accident, hydrogen explosion induced to a reduced public confidence in nuclear safety [6]. Other, hydrogen distribution, combustion, and mitigation studies have been applied in nuclear power plant models. Royl et al. [7] have made a hydrogen risk analysis during severe nuclear accident using CFD codes, to obtain localized detailed information and supplement the results of lumped parameter codes, which focus on global or average effect.

Oxidation of the cladding, rods, and other components in the core constructed in zirconium base alloy by steam is a critical issue in LWR accident producing severe core damage. The oxygen consumed by the zirconium is supplied by the upflow of steam from the water pool below the uncovered core, supplemented in the case of PWR by gas recirculation from the cooler outer regions of the core to hotter zones. In BWR, the gas recirculation process is prevented, since each fuel assembly is housed in its own channel box [8].

Fuel rod cladding oxidation is then one of the key phenomena influencing the core behaviour under high-temperature accident conditions. The chemical reaction of oxidation is very large exothermic, which determines the hydrogen rate generation and the cladding brittleness and degradation [9]. As the cladding material construction, Zircaloy-4 (Zry-4), oxidation process has been extensively studied for decades, whether in gaseous medium with air or oxygen or steam [14, 810]. Those studies showed that Zry-4 oxidation by air has similarities with oxidation in steam due to the common reaction partner oxygen, and also important differences. The exothermal heat released during air oxidation is around 1.8 times higher than steam, which causes a higher rise rate temperature [4].

Most of the works available in the literature dealing with the oxidation of Zry-4, agree that it is carried out in two temperature ranges. However experimental database at temperatures below 1800 K is very large and the oxidation kinetics in steam medium is well defined and their applicability. Variables predicted by them, as mass gain of oxygen and oxide growth, vary in the range of 20 to 30%. It is the same with which leads to higher temperatures of 1800 K, which is much more important uncertainty [1]. Since then, there have been several works to supplement the database at elevated temperatures. Schanz et al. [1] conducted a series of experiments to reduce uncertainty in the results and, lately, Shi and Cao [10] reports their work in compiling experimental data at elevated temperatures. Thus, in these works, there is reliable information that allows us to know the rate of hydrogen can be produced in the reactor core derived from the oxidation of Zry-4. However, there is still little information regarding the hydrogen generation rate as a function of the temperature profile in the fuel rods and the steam mass fraction. These parameters are very important to foresee the magnitude of accumulated hydrogen and the risk of this accumulation during a severe accident.

In this paper, we present the numerical analysis of heat and mass transfer with hydrogen generation in an array of fuel cylinder bars, each coated with a cladding and a current of water vapor flowing outside the cylinders. This analysis considers the fuel element without mitigation effects. The method applied consists in a representative periodic unit cell where the initial and boundary-value problems for heat and mass transfer (water reduction) were solved. The unit cell consists in a fuel element and steam water which allow analyzing the problem in small scale length. In this approximation, the fuel element and steam are assumed continuous, that is, the hypothesis for the heat flux in the interphase fuel-steam: where is the heat transfer coefficient, is the wall temperature and is the steam temperature far from the well, is not applied in this work. The paper is organized as follows: In Section 2, we provide an overview of the transport and reaction phenomena taking place in the system; in Section 3, we present the mathematical model that was numerically solved; the results are discussed in Section 4; finally in Section 5, we present the corresponding conclusions of the work.

2. Overview of Transport and Reaction Phenomena

Steam reduction occurs at the solid surface of the cladding and roads. Accordingly with Beuzet et al. [4], the oxidation process of the cladding and rods is divided into two regimes called pre- and postbreakaway regimes. During the prebreakaway regime, oxide growth is controlled by oxygen diffusion inside the oxide layer. Diffusion leads to equilibrium of the chemical species concentration. Assuming that hydrogen is present as a diluted solute in the steam flow allows applying Fick’s law as a constitutive equation for the diffusive mass flux. Then, the reaction can be divided in two steps: Then, the O/Zr ratio at the solid surface is a function of the steam-H2 ratio in the adjacent gas. Oxygen at the solid surface moves through the oxide scale by solid-state diffusion. At the oxide-metal interface, oxygen from steam reacts with the Zr to produce a substoichiometric oxide , that is, where the O/Zr ratio in the metal at the interface is the final solubility of oxygen in αZr to βZr. During this regime, oxide growth is controlled by oxygen diffusion inside the oxide layer. This phenomenon also follows Fick’s law so that the diffusion rate is proportional to the concentration gradient.

Then, the oxide layer loses its protectiveness for a critical thickness, whose value depends on the initial state of the cladding, the composition of the atmosphere, and temperature [3]. According to Duriez et al. [11], there are breaks at the surface layers of the oxide, due to the accumulation of stress and change of solid phase of Zr. In other words, there is an interrelationship between the change in density of the phase shift (αZr to βZr) and the difference of thermal expansion for the growth of the oxide layer and metal [4, 12]. Additionally, zirconia thus obtained is not pure, derived from the composition of Zry-4 and the rapid transition in the conditions that occurs during a severe accident: the crystalline phase changes of the oxide formed, monoclinic to tetragonal between 1133 K and 1473 K, and tetragonal to cubic between 1773 K and 2723 K [13]. Then, the oxide layer loses its protectiveness as a function of the cladding initial state (if an oxide layer is already exists), atmosphere composition, and temperature. Then, the stoichiometric oxidation reaction is The chemical equilibrium of this reaction shows a mechanism in which there exists a first-phase component of ZrO(OH)2 near to the ambient temperature.

At breakaway transition, the oxide rate increase after reached a minimum, and then inside the crack formed during this step, there is a fast oxidation inside the cladding, which contributed to form porosities [11]. The oxidizable surface is then increased for those cracks and also for the cladding distortion and so called ballooning. This last step is known as the postbreakaway regime [4].

In Figure 1, we show the computed chemical equilibrium data in the presence of ZrO2. We notice that the data is only stable at high temperatures. Hydrogen is produced at low temperatures, although with a very slow kinetics.

The heat of the reaction given by (1) is expressed as a function of temperature as: The transition solid phase enthalpy is  kJ/mol at 1445 K. These expressions for the enthalpy are used in the following section to take into account the heat generated at the solid-fluid interface.

The reaction of Zry with steam at elevated temperatures involves the growth of discrete layers of oxides and oxygen rich phase from parent -phase [1]. The main experimental data obtained in laboratory experiments includes oxide scale growth (), which accordingly with the chemical equilibrium is ZrO(OH)2 and total oxygen mass gain (). These quantities can be defined through where is the thickness of the layers, is the mass of oxygen absorbed by Zry per unit area, and is the reaction time.

Parabolic correlation is given as temperature-dependent Arrhenius-type function with constant activation energy and a factor : There are several classical parabolic correlations used to calculate the reaction rate in cladding oxidation at elevated temperatures. The kinetics evaluation of the growth of ZrO2 and layers was evaluated by several authors, nevertheless Urbanic and Heidrick [14] were the first to identify a discontinuity in the Arrhenius plot of the rate coefficients for mass increase and ZrO2 scale grows at the beginning of tetragonal to cubic transition zirconia. More recently, Shantz et al. [1] made a kinetics evaluation of ZrO2 and in the temperature range 1273–1773 K by calculated approximation of oxygen uptake comparing with some measurements. Other authors have also studied the kinetics of this reaction; nevertheless, we only consider those formulated by Leistikow and Schanz [15], Prater et al. [16], and Volchek [17].

The kinetics of the oxidation reaction is divided in three parts, in each one, a phase transition of the zirconium and zirconia from to and monoclinic to tetragonal and to cubic, respectively is involved. Leistikow and Schanz correlations were obtained from experimental sets, over the temperature range 973 to 1873 K. Weight gain and many specimens for growth of ZrO2 and were measured gravimetrically. These correlations correspond to a temperature range of 973 K to 1873 K: where is the parabolic coefficient of total oxygen mass gain (kg/m2s0.5), is the coefficient of oxide scale growth (m/s0.5).

According to Veshchunov et al. [18], the oxide layer consists of two sublayers in the temperature range of 1800 K to 2650 K: the tetragonal phase outside and the cubic phase inside. Later, the best-fit correlations were made at high temperature  K and verified the applicability to different of temperature transient based in its experiments [9] and also reported by Shi and Cao [10]. Then, the next equations can be applied to compute the kinetics parameter for the superficial oxidation that is predominant in these phenomena: Urbanic and Heidrick [14] correlation for mass gain rate have two turning points: at  K and  K, because the authors consider the tetragonal to cubic transition temperature as  K. Meanwhile, Fichot et al. [9] consider that the transition temperature is lower and at the beginning of the tetragonal to cubic transition  K oxidation behavior runs gently and rates increase smoothly. The oxidation rates calculated by Prater and Courtright [16] and Urbanic and Heidrick [14] agree with the best-fitted values and the errors are very small at  K. Nevertheless, the values calculated by Prater et al. [16] are better than other two correlations at  K [10]. For this reason, in this work, we consider the temperature ranges as  K; ;   K.

The oxygen flux available at the cladding surface strongly depends on the steam mass fraction in the bulk and properties of neighboring boundary layer: where is the mass transfer coefficient, and are steam mole fraction in the bulk and on the solid surface, respectively.

Assuming identity between Sherwood (Sh) and Nusselt (Nu) numbers, at outer cladding surface can be estimated as where is the diffusion coefficient in boundary layer, and is the equivalent hydraulic diameter of fluid channels associated with a single rod.

The maximum steam flux is obtained with , which means that all steam is consumed at the surface. Then, can be rewritten as According with Olander [8], for an axial steaming velocity of  m/s, ideal gas densities corresponding to typical total pressure and temperature in the event, an  m, the Reynolds number is lower than 2000. For this regime, the Nusselt number is 3.66. Then, it is possible to assume, as made by Olander, that is the binary diffusivity of steam in H2O(g)-H2 mixture. The expressions for these coefficients are [1] Densities (international units) used in this work of pure materials are [19] The heat capacity is computed by means of the coefficients of the polynomial formula [20]: , the coefficients for the materials used are presented in Table 1.

3. Mathematical Model in Reactor Fuel Elements

The system under consideration is depicted in Figure 2, which consists of an array of fuel cylinder bars, each coated with a cladding and a current of steam flowing outside the cylinders. In the same figure, we also show a representative periodic unit cell where the initial and boundary-value problems for heat and mass transfer are to be solved.

In this way, we can identify four regions of interest, where the problem under study has a characteristic length of the order of , that is, 16.2 mm, in contrast with the fuel assembly that is ten times larger as illustrated in Figure 3(b) (large scale). Accordingly, the unit cell regarding the core of the BWR is hundreds of times smaller (Figure 3(a)). Then, we can say that the study presented in this work takes place in a small length scale (Figure 3(c)).

Let us identify four regions in the unit cell, representative of the fuel element and coolant:(i) Region I: Fuel, ,(ii) Region II: Gap, ,(iii) Region III: Cladding, ,(iv) Region IV: Fluid (coolant), .

Since and continuity conditions of heat transfer can be assumed; we may neglect Region II overall. In this way, the differential equations for heat and mass transfer to be solved taking into account only three regions.

3.1. Heat Transfer Process

The fuel element temperature distribution was obtained considering each radial node at each of the twelve hydraulic axial nodes in the core. The fuel heat transfer formulation is based on the following fundamental assumptions: (i) axis-symmetric radial heat transfer, (ii) the heat conduction in the axial direction is negligible with respect to the heat conduction in the radial direction, (iii) the volumetric heat rate generation in the fuel is uniform in each radial node, (iv) storage of heat in the fuel cladding and gap is negligible, and (v) periodic boundary condition at entrances and exits of the unit cell.

Under these assumptions, the transient temperature distribution in the fuel element, initial and boundary conditions are given by each region:

Region I:

Region II: continuity conditions of heat transfer are assumed;

Region III:

Region IV:

Periodicity:

In these equations, the subscripts , and are used for denoting the fuel region, cladding region, and coolant region, respectively; is the temperature, is the density, is the heat capacity, is the conductivity, is the neutronic power per unit volume (which is calculated by standard law of the decay heat), and is the unit normal vector. We used the symbol to represent the volume-averaged variables, that is, is the average fuel temperature and is the average velocity of the coolant. Clearly, the averaging domains for and are not the same. Furthermore, (17) considers that the reaction heat due to chemical reaction between the Zr and coolant is represented by , which is given by (4) as a function of temperature. The periodic boundary conditions for unit cell which consider each regions is represented by (19), where and are constants and is a periodic function whose characteristic is zero mean [22]. Finally, at each section is assumed to be at different constant temperatures.

3.2. Mass Transfer Process

The mass transfer phenomena were discussed in Section 2, in this work, we consider that the hydrogen generated diffuses in the coolant by convection and diffusion. Then, the governing equation and boundary conditions for mass transfer are given by:

Periodicity:

where is the hydrogen concentration, is the diffusion coefficient, reaction rate coefficient, and and are constant. The hydrogen concentration at the initial time is assumed to be constants. Actually, for the developments that follow, the concentration is made dimensionless with this initial value.

The heat transfer processes are coupled with the mass transfer process through reaction rate and heat reaction at the interface, which are a functions of the temperature. Thus, in this work the simultaneous transfer of heat and mass is considered.

3.3. Decay Heat

Evaluation of the heat generated in a reactor after shutdown is important for determining cooling requirements under normal conditions and accident consequences. Reactor shutdown heat generation is the sum of heat produced from fission due to delayed neutron or photoneutron emissions, decay of fission products fertile materials, and activation products. The heat decay level used in this work is given by [23] where is the steady-state volumetric heat-generation rate, is time after shutdown, and is the time after reactor startup. Equation (22) is used in the heat transfer process in the fuel (14).

4. Numerical Experiments and Discussions

In order to solve the boundary-value problem presented in the previous section, we used the commercial finite-element solver Comsol multiphysics 4.2. The space and temporal meshing were adapted depending of the range of time of interest. Standard meshing and convergence analyses were performed in order to ensure the accuracy and exactness of the numerical results. The numerical solution was easily adapted to study each of the 12 nodes of the system. In the following paragraphs, we will present the numerical simulations for a particular node (node 3) that was arbitrarily chosen for illustration purposes. In order to gain some insight about the influence of the initial temperature distributions over the system transient performance, we computed four types of average temperatures, three corresponding to regions I, III, and IV and the fourth one for the whole unit cell. In the following paragraphs we discuss our results.

4.1. Average Temperatures

In Figure 4, we present the evolution of the average temperatures for the four types of averages mentioned above. We observe that the fuel temperature is practically unaffected by the heat transfer processes taking place in the other parts of the system. However, the average temperatures in the cladding and in the coolant exhibit more plausible oscillations in the temperature. It is interesting to notice that the average temperatures corresponding to the cladding and water regions fluctuate until reaching equilibrium at about 0.1 h. Furthermore, after the first hour is elapsed, the exothermic reaction and the heat source in the fuel make the average temperatures in the cladding and water regions to increase, and eventually the whole system will reach equilibrium. Overall, the whole system did not exhibit drastic temperature variations in the time range studied. This is to be expected since the maximum initial temperature differences are about 50 K. Let us also note that the average temperatures in the cladding and in the coolant regions tend to increase rapidly after about 10 hours. Certainly, if the heat source within the fuel section overcomes the heat generated by the chemical reaction at the solid-fluid interface, the steady state solution (not shown in the figure) should consist in a thermal equilibrium condition with respect to the fuel temperature.

In Figure 5, we compare the same average temperature transient profiles but imposing a larger initial heat flux driving force, that is, the initial temperatures in the cladding and in the coolant were taken to be half the values used in Figure 4. Now, we observe that the average temperatures corresponding to the cladding and the water exhibit a step-wise time evolution, reaching equilibrium at approximately the same time range as in the previous case. Interestingly, the values of the temperatures in these regions at the last computed times (corresponding to a 48 h simulation) approach the range of initial values shown in Figure 4. It is also worth noticing that the temperatures from the cladding and the coolant approach the one corresponding to the whole unit cell at a faster rate than in the previous case.

Another degree of freedom that we explored in our simulations was the fluid flow, which impacts the convective heat and mass transfer. Taking as reference the conditions used in Figure 4 and increasing the fluid velocity by a factor of 10, we obtained the results shown in Figure 6. In this case, the analysis had to be stopped at 3600 s due to problems in solution convergence generated by the large convective effects. Despite this drawback, we can observe that the time required for the first stabilization of temperatures is reached sooner than under the conditions used in Figure 4.

To gain a better insight of the transport phenomena taking place in the system, we provide in Figure 7 surface plots with the spatial temperature distribution for certain times. The conditions under which these figures were obtained correspond to those used in Figure 5. Here, we observe that the most dramatic changes of temperature take place near the boundaries of the system regions. For shorter times (i.e., at the first hour), the temperature variations are quite steep near the fuel and cladding boundary, whereas after ten hours these profiles become shallower. In addition, let us notice that the temperature profiles in the cooling steam region do not change very drastically with position but the same is not true for the time dependence. This means that the results from taking the spatial average of the temperature are justifiable since they provide, at the very least, a qualitative assessment of the physical phenomena taking place in the system.

4.2. Hydrogen Average Concentration

To finalize this section, we present in Figure 8 the average concentration profiles for the three cases considered above. In this context, Case A corresponds to the parameters used to obtain Figure 4, Case B for those in Figure 5, and Case C for those corresponding to Figure 6. The hydrogen concentration is displayed dimensionless with respect to its initial value. We observe that the initial temperature distributions do not play a role as relevant as that played by convection, since the characteristic times at which the concentration rapidly departs from its initial values are defaced in two orders of magnitude when the velocity is increased in only one order of magnitude. In other words, an increment in the velocity influences more drastically convective mass transfer than convective heat transfer. This is reasonable since heat transfer is mainly driven by the source in the fuel region, whereas mass transfer is driven by the interfacial reaction at the solid-fluid interface.

5. Conclusions

In this work, we carried out numerical simulations in a periodic unit cell that represents a fuel rod with its cladding and cooling steam of a BWR. Our main simulation variables were the initial heat distributions and the flow rate. We found that the temperature is more sensitive to changes on the initial distribution than over the flow rate. This follows from the fact that a larger temperature difference is established thus promoting heat transfer over time. Our main result is the prediction of the hydrogen generation that is carried with the cooling steam. In this case, we found that an increment in the flow rate resulted in a plausible decrement of the characteristic times at which the hydrogen generation exponentially increases. In addition, a change of the initial temperature distributions did not provide such drastic time variations.

These results can be used in further upscalings for decision making in terms of risk analysis. As shown in Figure 3, the real system under consideration is hierarchical in nature. This means that the results of this work can be used in conjunction with an appropriate upscaling technique in order to derive mathematical models at other levels of scale of the nuclear reactor. This type of modelling consisting of linking averaging and pointwise numerical simulations is currently known as hybrid upscaling and will be pursued in future works.