Heat and Mass Transfer during Hydrogen Generation in an Array of Fuel Bars of a BWR Using a Periodic Unit Cell

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.


Introduction
During a severe nuclear reactor accident, several works are focusing on core degradation by metal core components oxidation by air or steam [1][2][3][4].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 [1][2][3][4][8][9][10].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 Science and Technology of Nuclear Installations 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: q = h(T w − T ∞ ) where h is the heat transfer coefficient, T w is the wall temperature and T ∞ 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.

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-H 2 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 ZrO (2−x) , 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 ZrO 2 .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 ΔH t = 8.4 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 (k ox ), which accordingly with the chemical equilibrium is ZrO(OH) 2 and total oxygen mass gain (k T ).These quantities can be defined through where l ox is the thickness of the layers, W is the mass of oxygen absorbed by Zry per unit area, and t is the reaction time.Parabolic correlation k is given as temperature-dependent Arrhenius-type function with constant activation energy E a and a factor A: 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 ZrO 2 and α − Zr(O) 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 ZrO 2 scale grows at the beginning of tetragonal to cubic transition zirconia.More recently, Shantz et al. [1] made a kinetics evaluation of ZrO 2 and α − Zr(O) 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 ZrO 2 and α − Zr(O) were measured gravimetrically.These correlations correspond to a temperature range of 973 K to 1873 K: where K T is the parabolic coefficient of total oxygen mass gain (kg/m 2 s 0.5 ), K ox is the coefficient of oxide scale growth (m/s 0.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 T > 1800 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 K ox 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 T = 1853 K and T = 1873 K, because the authors consider the tetragonal to cubic transition temperature as T = 1853 K.Meanwhile, Fichot et al. [9] consider that the transition temperature is lower and at the beginning of the tetragonal to cubic transition T ≈ 1800 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 T > 1900 K. Nevertheless, the values calculated by Prater et al. [16] are better than other two correlations at T > 1900 K [10].For this reason, in this work, we consider the temperature ranges as T < 1800 K; 1800 < T < 1900; T > 1900 K.
The oxygen flux available at the cladding surface F gas strongly depends on the steam mass fraction in the bulk and properties of neighboring boundary layer: where k gas is the mass transfer coefficient, C bulk and C surf are steam mole fraction in the bulk and on the solid surface, respectively.
Assuming identity between Sherwood (Sh) and Nusselt (Nu) numbers, k gas at outer cladding surface can be estimated as where D g is the diffusion coefficient in boundary layer, and L h is the equivalent hydraulic diameter of fluid channels associated with a single rod.
The maximum steam flux is obtained with C surf = 0, which means that all steam is consumed at the surface.Then, F gas can be rewritten as According with Olander [8], for an axial steaming velocity of v steam ∼ = 0.1 m/s, ideal gas densities corresponding to typical total pressure and temperature in the event, an L h ∼ = 0.01 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 D g is the binary diffusivity of steam in H 2 O(g)-H 2 mixture.The expressions for these coefficients are [1] D H2 = 2.15 × 10 −5 T 2.33 P T , D bin = 6.262 × 10 −4 T 2.33 P T . ( Densities (international units) used in this work of pure materials are [19] βZr : ρ(T)=6434× 1 + 9.7 × 10 −6 ×T −3 , The heat capacity is computed by means of the coefficients of the polynomial formula [20]: Cp = A+BT +CT −2 +DT 2 , the coefficients for the materials used are presented in Table 1.

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 P, 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)).
Since r 2 /r 1 = 0.005321/0.005207= 1.021894 ≈ 1 and con- tinuity 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.

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: + T 0 , at entrances and exits of the unit cell. ( In these equations, the subscripts f , c, and f l are used for denoting the fuel region, cladding region, and coolant region, respectively; T is the temperature, ρ is the density, Cp is the heat capacity, K is the conductivity, q is the neutronic power per unit volume (which is calculated by standard law of the decay heat), and n is the unit normal vector.We used the symbol • to represent the volume-averaged variables, that is, T f f is the average fuel temperature and v is the average velocity of the coolant.Clearly, the averaging domains for T f f and v are not the same.Furthermore, (17) considers that the reaction heat due to chemical reaction between the Zr and coolant is represented by ΔH r , 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 h and T 0 are constants and (r) is a periodic function whose characteristic is zero mean [22].Finally, at t = 0 each section is assumed to be at different constant temperatures.

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: + c 0 , at entrances and exits of the unit cell, where c H2 is the hydrogen concentration, D is the diffusion coefficient, k reaction rate coefficient, and i and c 0 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.

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 q 0 is the steady-state volumetric heat-generation rate, (τ − τ s ) is time after shutdown, and τ is the time after reactor startup.Equation ( 22) is used in the heat transfer process in the fuel (14).

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.

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.
Science and Technology of Nuclear Installations 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 stepwise 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.

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.

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

Figure 3 :
Figure 3: Characteristic lengths of the system [21].(a) Nuclear reactor core; (b) Large scale (showing four fuel assemblies); (c) Small scale; σ-phase is the fuel element and γ-phase is the coolant.

Figure 4 :
Figure 4: Temporal evolution of the average temperatures in the system taking initial temperature distributions of 1115 K, 1126 K and 1169.5 K for the coolant, cladding, and fuel regions, respectively.The heat source term was fixed as 3403.569W/m 3 ; the average fluid velocity is 0.7045 m/s.

Figure 5 :
Figure 5: Temporal evolution of the average temperatures in the system taking initial temperature distributions of 557.5 K, 563 K and 1169.5 K for the coolant, cladding, and fuel regions, respectively.The heat source term was fixed as 3403.569W/m 3 ; the average fluid velocity is 0.7045 m/s.

Figure 6 :
Figure 6: Same case as Figure 4 but with the velocity increased 10 times.

Figure 7 :
Figure 7: Evolution of the temperature profiles in the unit cell for the conditions used in Figure 5.