Heat Transfer and Hydraulic Resistance in Nuclear Fuel Rods

Abstract

This chapter briefly describes the fundamental concepts of heat transfer and hydraulic resistance in water-cooled nuclear reactors, more specifically the nuclear fuel assemblies. There are two key areas in nuclear thermal-hydraulics, namely heat transfer and hydraulic resistance with and without phase change. Boiling and condensation play a critical role in water-cooled nuclear reactors as these are needed for the design, operation and safety analysis of nuclear reactors. The common models used in the nuclear industry are described in this chapter.

1. Introduction

The word thermal-hydraulics is composed of two Greek terms which indicate heat and water respectively. Therefore, thermal-hydraulics studies the behaviour of fluid(s) subjected to heat. In particular, nuclear thermal-hydraulics is more specific to the analysis of nuclear systems, such as nuclear power plants (NPPs). It is used mainly for the design of thermal systems, dimensioning of thermal or hydraulic components, system performance and safety analysis. Key phenomena associated with this discipline are boiling heat transfer, hydraulic resistance under single-phase and multi-phase conditions.

The field of nuclear thermal-hydraulics, specifically for water-cooled reactors, improved significantly during the peak of the nuclear industry. This resulted in the production and publication of a large body of knowledge, such as books, articles, conference proceedings. It would be impossible, and impractical, to replicate all this information in a single chapter. For that reason, this chapter entitled “Heat transfer and hydraulic resistance in nuclear fuel rods” presents only the fundamentals of nuclear thermal-hydraulics, with a focus on CANDU (CANada Deuterium Uranium) type reactors. For that reason, the model developments and equations are omitted. The reader is instead led to make use of the references cited in this chapter for further details.

2. Reactor power

Important input needed for thermal-hydraulics analysis is the core power distribution, that is, how the power is distributed within a reactor’s core. This analysis is usually determined by carrying out reactor physics calculations. In CANDU reactors, the fuel management differs from pressure vessel water-cooled reactors, because it uses natural uranium, which requires a different fuel management process. On-power refuelling is used to replace fuel as needed, to keep the reactor operating as designed. Furthermore, CANDU reactors use independent fuel channels to contain their nuclear fuel, and consequently, each fuel channel produces different powers.

To start the reactor power analysis, the link between the neutron flux with power is represented by the following equation:

$$q^{\prime \prime \prime }\left(r\right)=E_d{\int }_0^{\infty }{\Sigma }_f\left(E\right)\Theta \left(r,E\right)\mathit{dE}$$E1

Where E_d is the energy deposited locally in the fuel per fission, ${\mathrm{\Sigma }}_f$ is the macroscopic fission cross-section of the fuel, and $\mathrm{\Theta }$ is the neutron flux. This equation essentially means that power generated per unit of volume is the neutron reaction rate times the energy deposited per fission. For a cylindrical homogeneous reactor without a reflector, the neutron flux can be approximated [1].

$$\Theta \left(r,z\right)=A{J}_o\left[\frac{2.405r}{R_e}\right]\mathit{cos}\left[\frac{\mathit{\pi z}}{H_e}\right]$$E2

Where:

$$A=\frac{3.63P}{V{E}_R{\mathrm{\Sigma }}_f}$$E3

$R_e$ and $H_e$ are the extrapolated lengths of the radius and height of the cylinder, and J_o is the Bessel function of order zero, $E_R$ is the recoverable energy in joules per fission and $V$ the volume of the cylinder.

After some mathematical manipulations, the rate of heat production per unit of volume of a fuel rod is [1, 2]:

$$q^{\prime \prime \prime }\left(r,z\right)=q_{\mathit{\text{center}}}^{\prime \prime \prime }{J}_o\left[\frac{2.405r}{R_e}\right]\mathit{cos}\left[\frac{\mathit{\pi z}}{H_e}\right]$$E4

The solution of this equation is depicted in Figure 1. This figure shows the axial $\left(z\right)$ and radial $\left(r\right)$ power distribution.

Figure 1.

Representation of an axial and radial power distribution.

In reality, the power distribution is an implicit calculation between neutron physics and thermal-hydraulics. In order to compute the neutron flux and power, it is necessary to know the coolant and fuel temperature distributions, which are unknown until thermal-hydraulics analyses are performed. However, the latter analyses require power and power distribution, and this results in an implicit calculation. Usually, a coupled and iterative process is used to get a solution. The reader is referred to reference [2] for additional information.

From Eq. (4), it can be stated that the power distributions are non-uniform, in fact, these have a cosine-like profile. The fact that the radial (core) power distribution is non-uniform has important effects on the dimensioning of equipment, systems and/or components of a nuclear core. For example, in CANDU reactors some fuel channels produce more power than others. The ones near the centre usually produce more than those that are in the periphery. Furthermore, as the inlet and outlet coolant temperatures are fixed for all the fuel channels, it is important to ensure each fuel channel has received the right amount of coolant to remove the heat generated by the nuclear fuel. Indeed, fuel channels with lower power require less mass flow than fuel channels with higher power. In CANDU reactors, the channel mass flow is not a controlled parameter. Orifices are used in the fuel channel to provide the right amount of coolant in the fuel channels. However, under transient conditions or accidents, the channel flow varies with changes in the reactor and channel powers, and with the thermal-hydraulics conditions at the inlet and outlet headers.

As CANDU reactors are refuelled online, the insertion of fresh bundles in a channel, which usually replaces four bundles at a time-creates a local-power effect due to the inserted reactivity (see Figure 2). This is called ripple effect and can impact several neighbouring fuel channels. In addition, the continuously changing burn up distribution can result in power oscillations due to Xenon or regional overpower¹. These overpower limits are separate from and above the normal operating limits on channel and bundle powers.

Figure 2.

Schematic of a CANDU power distribution.

Another important characteristic of CANDU fuel channels is the End-Flux-Peaking (EFP). This is a phenomenon that affects the flux profile of a fuel bundle. It occurs in the end region separating two individual bundles. The geometry in the end regions consists of D₂O coolant, a Zircaloy endplate and end caps, and uranium dioxide fuel pellets. Due to the fact that the heavy-water-coolant and Zircaloy structural material have a much lower absorption cross-section for thermal neutrons than uranium dioxide, thermal neutrons tend not to be absorbed as much as they would be in the fuel. The lack of absorption of thermal neutrons in these regions leads to peaks in neutron flux around the end regions. The occurrence of EFP leads to higher fission rates, leading to more heat being produced and therefore overall higher temperatures within the pellets that are adjacent to these end regions. Higher temperatures lead to an increased risk of sheath strain, corrosion, and fuel centreline melting, all of which have a significant impact on the integrity of the sheath and fuel pellets during an overpower of loss-of-coolant-accident conditions.

2.1 Heat characteristics of a CANDU 6 nuclear reactor

In a nuclear reactor, the heat produced by the nuclear reaction is primarily removed through conduction and convection. For example, the heat produced in the nuclear fuel is transferred by conduction to the surface of the rod. Convection involves the transfer of heat by the movement of a fluid. Thus, the heat conducted to the surface of a fuel rod is transferred to the coolant by convection. In addition, under some specific accidents, where the fuel reaches high temperatures, radiation heat transfer plays an important role.

2.1.1 Heat conduction equation

Fortunately, as the reactor core, fuel channels and fuel elements have a long cylindrical shape, the heat conduction equation can be written in cylindrical coordinates, which can then be used as a base for the heat conduction analysis. The general heat conduction equation is:

$$\rho C_p\left(r,T\right)\frac{\partial T}{\partial t}=\nabla \bullet k\left(T\right)\nabla T+q^{\prime \prime \prime }$$E5

Where $C_p$ is the specific heat at constant pressure (which for incompressible materials is equal to the specific heat and constant volume $C_v$), $\rho$ is the density of the fuel, $k$ is the thermal conductivity, $T$ is the temperature and $q^{\prime \prime \prime }$ is the volumetric heat.

Note that the conductivity of the fuel is temperature dependent, as shown in Eqs. (6) and (7), which are formulations for the thermal conductivity of Uranium [3]:

$$k=\frac{100}{6.548+23.533T}+\frac{6400}{T^{5/2}}exp\left(\frac{-16.35}{T}\right)$$E6

Where $T$ is defined as $T/1000$, and $T$ is in Kelvin. The conductivity $k$ is the thermal conductivity for 95% dense UO₂ in W/m K.

Eq. (6) can be approximated using polynomials [3]:

$$\begin{array}{l}k=+12.57829-2.31100\times {10}^{-2}T+2.36675\times {10}^{-5}{T}^2-1.30812\times {10}^{-8}{T}^3\\ +3.63730\times {10}^{-12}{T}^4-3.90508\times {10}^{-16}{T}^5\end{array}$$E7

To better understand the analysis of the heat conduction equation, the following example is that of a classical heat transfer problem used in nuclear reactor analyses, which is the analysis of a single fuel element (see Figure 3 for details). The objective is to solve the heat conduction equation, in order to determine the temperature profile within the fuel element. Both the maximum centreline temperature and cladding temperature are needed to ensure the integrity of the fuel and cladding is not compromised. The average temperature of the fuel is also needed for nuclear physics calculations. For example, assuming a predominant radial heat transfer under steady-state conditions, which is common in a nuclear fuel rod because the cylinder length is several times larger than the diameter, Eq. (5) is recast as:

Figure 3.

Cross-section geometry of a fuel element.

$$\frac{1}{r}\frac{\partial }{\partial r}\left(rk\frac{\partial T}{\partial r}\right)+q^{\prime \prime \prime }=0$$E8

Two boundary conditions are needed to solve this equation:

1. A boundary condition specifying that the heat flux at the centre of the fuel element is zero, as it is axisymmetric with respect to the axial direction,

2. The description of the heat transfer between the bulk temperature of the fluid $\left(T_b\right)$ and the outer cladding (wall) temperature of the element $\left(T_w\right)$ using Newton’s law of cooling.

$${\left.\frac{\partial T}{\partial r}\right]}_{r_0=0}=0\mathit{\text{for}}r=0$$E9

$$\begin{array}{ccc}{\left.q^{"}=-k\frac{\partial T}{\partial r}\right]}_{r=w}=\stackrel{\sim }{h}\left(T_w-T_b\right)& \mathit{\text{for}}& r=w\end{array}$$E10

Fortunately, this equation is not difficult to solve. For simplified cases, it is possible to find an analytical solution. As the problem is more elaborate, the equation is frequently solved by using numerical methods. In reality, the problem is much more challenging, as the materials—such as the nuclear fuel—undergo several nuclear and chemical reactions that change their properties over time. In addition, the mechanical properties are also affected by the conditions to which the materials are exposed during their lifetime. Changes in geometry are also possible, such as element bowing or ballooning. In other cases, oxide deposits (crud) in the cladding need to be considered.

The nuclear industry uses specific codes to assess in detail the heat transfer from the fuel element to the coolant. These codes must take into account the change of properties of the materials over time as the fuel is consumed and fission products are produced inside the fuel. Chemical reactions and mechanical stresses that can change the geometry of the system should also be considered for a complete analysis of the fuel element.

2.2 Heat transfer to coolants

Thus far, the value of the heat transfer coefficient $\left(\stackrel{\sim }{h}\right)$ has not been discussed in detail. Calculating and/or accurately estimating this coefficient is one of the most important and more difficult parts in thermal-hydraulics analyses. The reason for that is that this coefficient depends on multiple variables, such as geometry, flow conditions, flow regime and coolant properties.

Note that the second boundary condition needs the value of $\stackrel{\sim }{h}$. This boundary condition links the conduction with the convection heat transfer mechanisms. In other words, the heat from a solid is removed by a fluid (liquid or gas). This is described according to Newton’s law of cooling [1]:

$$q^{\prime \prime }=\stackrel{\sim }{h}\left(T_b-T_w\right)$$E11

Where $q^{”}$ is the heat flux, $T_w$ is the temperature of the surface of the solid (in the case of reactor fuel, this is the outer temperature of the cladding or wall) and $T_b$ is the temperature of the fluid.

The numerical value of the heat transfer coefficient, $\stackrel{\sim }{h}$ is a function of several variables such as physical properties of the fluid, mass flow rate, geometry and orientation of the system.

The solution of Newton’s law of cooling for single-phase conditions has been commonly solved using the method of similarities, or dimensional analysis [4]. Three dimensionless numbers are used to correlate the heat transfer coefficient, namely the Reynolds $\left(\mathbf{Re}\right)$, the Nusselt $\left(\mathbf{Nu}\right)$ and Prandtl $\left(\mathbf{Pr}\right)$ numbers

$$\mathbf{Re}\equiv \frac{G{D}_e}{\mu }$$E12

$$\mathbf{Pr}\equiv \frac{C_p\mu }{k}$$E13

$$\mathbf{Nu}\equiv \frac{\stackrel{\sim }{h}D}{k}$$E14

The solution takes the form of:

$$\mathbf{Nu}=a{\mathbf{Re}}^b{\mathbf{Pr}}^c$$E15

Where $a$, $b$ and $c$ are coefficients, which are usually obtained experimentally.

For example, one of the most common correlations is the Dittus-Boelter correlation [5]:

$$\frac{D\stackrel{\sim }{h}}{k}=\mathbf{Nu}=0.023{\mathbf{Re}}^{0.8}{\mathbf{Pr}}^{0.4}$$E16

Another common correlation is the Sieder and Tate correlation [6]

$$\frac{D\stackrel{\sim }{h}}{k}=\mathbf{Nu}=0.023{\mathbf{Re}}^{0.8}{\mathbf{Pr}}^{\frac{1}{3}}{\left(\frac{{\mu }_b}{{\mu }_w}\right)}^{0.14}$$E17

Both correlations are applicable under single-phase conditions. However, the Sieder and Tate correlation [6] (Eq. (17)) uses an additional term to take into account the effect of the wall temperature, by also taking into account the viscosity evaluated at the wall temperature.

The boiling process is a clear example of a thermal-hydraulics analysis. It is evident that this process involves fluid dynamics and heat transfer. Fluid dynamics are needed to determine the momentum and enthalpy characteristics of the flow, such as velocity, pressure drop and void profiles. Heat transfer is also required to compute the temperature distributions of the system and heat transfer rates. The calculations are highly inter-related as one feeds back to the other. For example, the velocity of the flow affects the heat transfer rate, which in turn affects void generation (boiling) affecting the hydraulic resistance, which in turn can change the flow regime, and so forth, as these disciplines are all highly coupled. The next Section 3 briefly explains the forced flow boiling process.

3. Boiling heat transfer

As CANDU reactors are water-cooled, it is important to briefly describe the boiling process. This process is usually divided into two types: (1) pool boiling, in which there is no liquid mass flow entering or leaving the system, similar to that of an actual pool, and (2) the second deals with systems where the liquid mass flow enters and leaves the system. In a CANDU power plant, these two types of analysis are relevant. However, in this chapter, we focus solely on the analysis of a CANDU fuel channel. In a CANDU fuel channel, the coolant mass flow is forced as it is driven by the main heat transport pumps. As the liquid coolant moves through the fuel channel, the liquid absorbs the heat produced by the fuel bundles/elements from the nuclear reaction. As the coolant moves along the channel, it develops velocity and temperature distributions. These distributions dictate the process of boiling.

3.1 Boiling curve (forced flow)

Understanding the phenomenon of heat transfer under boiling conditions, or heat transfer with phase change, is critical in thermal-hydraulics analyses. A detailed review of the boiling process is provided in several nuclear engineering and convection heat transfer textbooks [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18].

The boiling process can be explained by assuming a system where a heat flux from a heated fuel rod is transferred to a liquid coolant, such as water. The results are measured as a function of the temperature of the surface of the rod for a given system pressure and flow rate. The results of this experiment are shown in Figure 4.

Figure 4.

Boiling curve.

The heat flux increases slowly as the rod temperature is increased at low values. In this temperature range, between points A and B, heat is transferred to the coolant by convention without phase change. The heat transfer coefficient is usually determined by empirical correlations, such as Eqs. (16) and (17).

As the surface temperature of the fuel is increased further, a point is eventually reached where bubbles or vapour form at various imperfections on the surface of the fuel rod. This occurs at about point B (in Figure 4). This type of boiling is called nucleate boiling. As the bubbles are formed, they are entrained from the rods into the bulk of the coolant as a result of the turbulent movement of the fluid. However, if the bulk temperature of the coolant is lower than its saturation temperature, the vapour condenses in the liquid disappearing from the coolant. Thus, there is no net production of steam under these circumstances. This boiling process is known as subcooled nucleate boiling or local boiling. If and when the bulk temperature of the coolant reaches its operation temperature, the bubbles remain within the coolant stream, there begins to be a net production of steam and the system is undergoing a saturated nucleate boiling or bulk boiling.

In any event, with the onset of nucleate boiling, the heat moves into the liquid. At every temperature in this region between B and C, heat transfer is more efficient than ordinary convection. There are two reasons for this. First, heat is removed from the rods both as heat of vapourization and sensible heat. Second, the motions of the bubbles lead to rapid mixing of the fluid. The rapid increase in the heat flux with temperature is explained by the fact that the density of the bubbles forming at and departing from the rod surface increases rapidly with surface temperature.

With increased vapourization in a coolant channel, the heated surface becomes intermittently exposed to patches of vapour. Since the heat transfer coefficient decreases when the surface is covered with vapour, the wall temperature rises correspondingly. Hence, the wall temperature may become unstable as the surface is alternately covered with vapour or liquid, but then rise after the wall liquid is completely vapourized. Such behaviour, characterised by a marked temperature rise of the heated surface during boiling, as a result of a change in the heat-transfer mechanism, is called boiling crisis. There are two types of boiling crisis, departure of nucleate boiling or DNB and dryout.

With the onset of the boiling crisis, the heat flux into the coolant begins to drop. This is due to the fact that over the regions of the rods covered by vapour film, the heat is forced to pass through the vapour into the coolant by conduction and radiation, both are relatively inefficient heat transfer mechanisms. The heat transfer continues to drop more or less erratically (dotted line) with increasing fuel temperatures as the total area of the film covering the fuel increases. In this region, the system is said to be experiencing partial film boiling.

Eventually, when the rod surface temperature is high enough, the vapour film covers the entire rod and the heat flux to the coolant falls to a minimum value (point D). Beyond this point, any increase in temperature leads to an increase in the heat flux simply because heat transfer through the film, although a poor and inefficient process, nevertheless increases with the temperature difference across the film. The system is said to be undergoing full film boiling.

The critical heat flux, at which burnout is expected to occur, is an important design consideration in water-cooled reactors. The knowledge of burnout conditions is important not only for the design of a fuel bundle and its maximum operating conditions at nominal power but also under upset conditions, such as might arise from loss of coolant flow conditions due to power excursions.

4. Fluid dynamics

As previously mentioned in Section 2, the energy (heat) generated in the nuclear reactor core is removed by the coolant mass flow. The thermal analysis of a nuclear power plant requires the knowledge of the pressure and velocity distribution of the coolant at different locations of the reactor, and the transfer of heat between fuel rods to the coolant, or between the primary and secondary sides.

There is a considerable amount of literature on fluid dynamics. For that reason, only the equations needed to explain a relevant topic in a nuclear power plant are hereafter provided. The solution or procedures to solve these equations are beyond the scope of this chapter.

4.1 Single-phase flow

Various forms of the single-phase fundamental equations exist. This chapter follows the format used by Delhaye et al. [15]. In this subsection, we present the fundamental conservation equations. These equations can be written in a general form as:

$${\int }_{V\left(t\right)}\frac{\partial \phi }{\partial t}\mathit{dV}-{\int }_{V\left(t\right)}\left(\overrightarrow{\nabla }\bullet \phi \overrightarrow{v}\right)\mathit{dV}+{\int }_{V\left(t\right)}\overrightarrow{\nabla }\bullet \overline{\overline{J}}\mathit{dV}-{\int }_{V\left(t\right)}{S}_g\mathit{dV}=0$$E18

Where $\phi$ is a vector containing the variables which are conserved, for example, mass, momentum and energy, per unit of volume. $\overline{\overline{J}}$ is the flow of property per unit of area and time across the surface that bounds the material volume $V\left(t\right)$, and $S_g$ is the generation of the property $\mathrm{\phi }$ per unit of volume and time.

This equation can be rewritten in terms of a partial derivative [15]:

$$\frac{\partial \phi }{\partial t}+\overrightarrow{\nabla }\bullet \phi \overrightarrow{v}+\overrightarrow{\nabla }\bullet \overline{\overline{J}}=S_g$$E19

The first term represents the time rate of change of the property $\phi$, per unit of volume, the second term is the rate of convection per unit of volume, the third term is the surface flux and the fourth term is the volume source.

Therefore the local continuity, momentum and energy equations can be summarised using the following vectors. The first element of the vector $\phi$ is the density (mass per unit of volume), the second row is the momentum and the third element is the energy.

$$\phi =\left[\begin{array}{c}\rho \\ \rho \overrightarrow{v}\\ \rho \left(v+\frac{1}{2}\overrightarrow{v}\overrightarrow{v}\right)\end{array}\right];J=\left[\begin{array}{c}0\\ P\overline{\overline{I}}-\overline{\overline{\sigma }}\\ {\overrightarrow{q}}^{\prime }-\overline{\overline{T}}\overrightarrow{v}\end{array}\right];S=\left[\begin{array}{c}0\\ \rho \overrightarrow{g}\\ \mathit{\rho g}\overrightarrow{v}+{\dot{Q}}_g\end{array}\right]$$E20

The stress tensor $-\overline{\overline{T}}$ represents the normal and shear stresses acting on the surface $\overline{\overline{J}}=-\overline{\overline{T}}=-\left(-p\overline{\overline{I}}+\overline{\overline{\sigma }}\right)$.

These are the instantaneous conservation equations for a single-phase flow. The total number of unknowns is six: velocity vectors (in three directions, thus three unknowns), pressure, temperature and density. An equation of state is also used to close the system of equations.

However, if the flow receives enough heat, it can undergo a phase change. For example, in a CANDU fuel channel, the liquid coolant can change to vapour as it removes the heat generated in the nuclear reactor. In this situation, the flow is said to be a two-phase flow.

4.2 Two-phase flow

There are several definitions of two-phase flow in literature. Yet, one of the most practical is provided by Shire (as cited by Butterwort [8]).

“Two-phase flow is a term covering the interaction flow of two phases (gas, liquid or solid) where the interface between the phases is influenced by their motion. The proviso concerning the interface is inserted to distinguish between problems which are usually considered as two phases and those which are normally accepted as a single phase.”

There are more definitions of two-phase systems that explain the difference between multi-component and multi-field systems [14]. This chapter is concerned only with water and vapour two-phase flows.

An important area of nuclear thermal-hydraulics is the two-phase system analysis. Furthermore, under certain specific conditions, such as in a loss of coolant accident (LOCA), additional components such as non-condensable gases can be part of the system, thus adding one component more to the liquid-vapour mixture. In this case, the system is said to be multi-phase and multi-component. The system can be further classified, namely, if the system presents different structural formations, such as liquid droplets or mist, it is also classified as multi-field.

The general local conservation equation takes the following form [15]:

$$\sum _{k=1}^2{\int }_{V\left(t\right)}\left(\frac{\partial {\phi }_k}{\partial t}-\left(\overrightarrow{\nabla }\bullet {\phi }_k{\overrightarrow{v}}_k\right)+\overrightarrow{\nabla }\bullet {\overline{\overline{J}}}_k-S_{\mathit{gk}}\right)\mathit{dV}=\underset{A_i\left(t\right)}{\int }\sum _{k=1}^2\left(\overrightarrow{u_k}\bullet \phi \left(\overrightarrow{v_k}-\overrightarrow{v_i}\right)+{\overrightarrow{u}}_k\bullet {\stackrel{̿}{j}}_k\right)\mathit{dA}$$E21

$$\phi =\left[\begin{array}{c}\rho \\ \rho \overrightarrow{v}\\ \rho \left(v+\frac{1}{2}\overrightarrow{v}\overrightarrow{v}\right)\end{array}\right];J=\left[\begin{array}{c}0\\ P\overline{\overline{I}}-\overline{\overline{\sigma }}\\ {\overrightarrow{q}}^{\prime }-\overline{\overline{T}}\overrightarrow{v}\end{array}\right];S=\left[\begin{array}{c}0\\ \rho \overrightarrow{g}\\ \mathit{\rho g}\overrightarrow{v}+{\dot{Q}}_g\end{array}\right]$$E22

Where the subindex $k$ represents a phase, gas $\left(g\right)$ or liquid $\left(l\right)$.

Note that the first term of Eq. (21) is similar to the single-phase general conservation equation (Eq. (19)). However, the introduction of the summation means that there are two terms, one for the liquid phase and another for the gas phase. Similar to Eq. (19), each of these terms represents the conservation of a property $\phi$. The right side term of Eq. (21) is the interface condition, that is, the properties that are transported from one phase to another via an interface.

Note that there are two velocities for each phase. Under two-phase flow conditions, the lighter phase tends to travel faster than the heavier phase. Similarly, the two-phase mixture can also be in thermodynamic non-equilibrium, that is, each phase can have independent pressure and temperature.

Solving these systems of equations is a very complex task, and is usually prohibitive. Fortunately, this approach is not needed, and the prediction of averaged quantities is usually sufficient.

Averaging procedures have been proposed by numerous researchers, such as Delhaye et al. [15] and Ishii and Hibiki [16].

4.3 Practical two-phase models

As stated in the previous sections, two-phase flow systems are defined by the presence of an interface separating the phases. A more complete description of a two-phase flow system can be obtained using the two-fluid model, which requires solving the conservation equations of mass, momentum and energy for each phase. This system of equations is complex due to the non-linearity nature of fluid dynamics.

To alleviate this complexity, some assumptions have been made to decrease the number of variables, and therefore the number of equations. However, the assumptions made should be in accordance with the nature of the problem. Fortunately, in a nuclear power plant, as in many other plants, the fluid is flowing through pipes where the fluid is essentially axially predominant, resulting in a one-dimensional (1-D) equations. This significantly reduces the complexity of the equations, as these can be written assuming 1-D flow. Nevertheless, there are components or systems where 2-D or 3-D effects should be taken into account, such as the detailed thermal-hydraulics analysis of a fuel bundle.

Some of the most common two-phase models are described in detail in several two-phase flow textbooks [2, 7, 8, 10, 11, 12, 13, 14, 15, 16, 17]. For that reason, only the conservation equations are herein presented. The most common models are:

1. Homogeneous equilibrium model;

2. Simplified two-fluid model for separated flow;

3. Drift flux model.

4.3.1 Homogeneous equilibrium model

The Homogeneous Equilibrium Model (HEM) is the simplest model. It assumes thermo-equilibrium between phases, as well as equal velocities between them. That is, the conservation equations used for single-phase flow are also applicable in this model, however, the density (${\rho }_m)$ and enthalpy ($h_m)$ are mixture properties.

This model is applicable if one of the phases is finely dispersed, and in which the momentum and energy transfers are sufficiently rapid for the average velocities and the average temperatures of the two phases to be equal.

The conservation equations for the HEM are:

1. Mass conservation equation

   $$\frac{\partial {\rho }_m}{\partial t}+\frac{\partial }{\partial z}{\rho }_{m}v=0$$E23

2. Momentum conservation equation

   $$\frac{\partial }{\partial t}{\rho }_{m}v+\frac{\partial }{\partial z}\left({\rho }_m{v}^2\right)+\frac{\partial P}{\partial z}+{\rho }_{m}g=-\frac{4{\tau }_w}{D}$$E24

3. Total energy conservation equation (enthalpy)

   $$\frac{\partial }{\partial t}{\rho }_m\left[h_m+\frac{G^2}{2{\rho }_m^2}+\mathit{gz}\right]-\frac{\partial P}{\partial t}+\frac{\partial }{\partial z}\left[G\left(h_m+\frac{G^2}{2{\rho }_m^2}+\mathit{gz}\right)\right]={q^{\prime \prime \prime }}_w$$E25

4.3.2 Simplified two-fluid model: Separated flow

This is a simple case of a two-fluid model. In this model, it is assumed that the gas and liquid phases are separate, such as in an annular flow as shown in Figure 5. The important assumptions for this model are:

• The space and time correlation coefficients are equal to 1.

• The pressure is constant over the cross-sectional area.

• The equation of state valid for local quantities applies to averaged quantities.

• Longitudinal conduction terms in each phase, as well as their derivatives, are negligible.

Figure 5.

Control volume for the separate flow model.

Applying these assumptions, the conservation equations for this model are:

1. Mass conservation:

   $$\frac{\partial }{\partial t}\mathit{A\alpha }{\rho }_g+\frac{\partial }{\partial z}\left(\mathit{A\alpha }{\rho }_g{v}_g\right)=\frac{1}{\mathit{dz}}\mathit{\delta w}$$E26
   $$\frac{\partial }{\partial t}A\left(1-\alpha \right){\rho }_l+\frac{\partial }{\partial z}\left(A\left(1-\alpha \right){\rho }_l{v}_l\right)=-\frac{1}{\mathit{dz}}\mathit{\delta w}$$E27

2. Momentum conservation equation:

   $$\frac{\partial }{\partial t}\mathit{A\alpha }{\rho }_g{v}_g+\frac{\partial }{\partial z}\mathit{A\alpha }{\rho }_g{v}_g^2+\mathit{A\alpha }\frac{\partial P}{\partial z}+{\mathit{A\alpha \rho }}_{g}g=\frac{\mathit{\delta v}}{\mathit{dz}}{v}_{g_i}-\frac{1}{\mathit{dz}}{A}_i\mathit{sin}\gamma {\tau }_g$$E28
   $$\begin{array}{l}\frac{\partial }{\partial t}A\left(1-\alpha \right){\rho }_l{v}_l+\frac{\partial }{\partial z}A\left(1-\alpha \right){\rho }_l{v}_l^2+A\left(1-\alpha \right)\frac{\partial P}{\partial z}+A\left(1-\alpha \right){\rho }_{l}g\\ =-\frac{\mathit{\delta v}}{\mathit{dz}}{v}_{l_i}+\frac{1}{\mathit{dz}}{A}_i\mathit{sin}\gamma {\tau }_g-\frac{1}{\mathit{dz}}{A}_w\mathit{sin}\theta {\tau }_w\end{array}$$E29
   $$\frac{\partial }{\partial t}\left[\alpha {\rho }_g{v}_g+\left(1-\alpha \right){\rho }_l{v}_l\right]+\frac{1}{A}\frac{\partial }{\partial z}A\left[\alpha {\rho }_g{v}_g^2+\left(1-\alpha \right){\rho }_l{v}_l^2\right]+\frac{\partial P}{\partial z}+{\rho }_{m}g=-\frac{A_w\mathit{sin}\theta {\tau }_w}{V}$$E30

3. Total energy conservation equation (enthalpy):

$$\begin{array}{c}\frac{\partial }{\partial t}\mathit{A\alpha }{\rho }_g{h}_g^{o-}\frac{\partial }{\partial t}\mathit{A\alpha P}+\frac{\partial }{\partial z}A\alpha {\rho }_g{h}_g^o{v}_g=-\frac{1}{\mathit{dz}}{\int }_{A_i}{m}_g{h}_g^o\mathit{dA}+\frac{1}{\mathit{dz}}{\int }_{A_i}{m}_g\frac{P}{{\rho }_g}\mathit{dA}-\frac{1}{\mathit{dz}}{\int }_{A_i}P{v}_g\left({\overrightarrow{u}}_g\bullet {\overrightarrow{u}}_z\right)\mathit{dA}\\ +\frac{1}{\mathit{dz}}{\int }_{A_i}{\overrightarrow{u}}_z·\left({\overrightarrow{u}}_g·\stackrel{̿}{{\sigma }_g}\right)v_g\mathit{dA}-\frac{1}{\mathit{dz}}{\int }_{A_i}\left({\overrightarrow{u}}_g\bullet {\overrightarrow{q}}_g^{\prime }\right)\mathit{dA}-\frac{1}{\mathit{dz}}{\int }_{A_{\mathit{gw}}}{\overrightarrow{u}}_{\mathit{gw}}·{\overrightarrow{q}}_l^{\prime }\mathit{dA}\end{array}$$E31

$$\begin{array}{c}\frac{\partial }{\partial t}A\left(1-\alpha \right){\rho }_l{h}_l^{o-}\frac{\partial }{\partial t}A\left(1-\alpha \right)P+\frac{\partial }{\partial z}A\left(1-\alpha \right){\rho }_l{h}_l^o{v}_l=-\frac{1}{\mathit{dz}}{\int }_{A_i}{m}_l{h}_l^o\mathit{dA}+\frac{1}{\mathit{dz}}{\int }_{A_i}{m}_l\frac{P}{{\rho }_l}\mathit{dA}\\ -\frac{1}{\mathit{dz}}{\int }_{A_i}P{v}_l\left({\overrightarrow{u}}_l\bullet {\overrightarrow{u}}_z\right)\mathit{dA}+\frac{1}{\mathit{dz}}{\int }_{A_i}{\overrightarrow{u}}_z·\left({\overrightarrow{u}}_l·\stackrel{̿}{{\sigma }_l}\right)v_l\mathit{dA}-\frac{1}{\mathit{dz}}{\int }_{A_i}\left({\overrightarrow{u}}_l\bullet {\overrightarrow{q}}_l^{\prime }\right)\mathit{dA}-\frac{1}{\mathit{dz}}{\int }_{A_w}{\overrightarrow{u}}_{\mathit{lw}}·{\overrightarrow{q}}_l^{\prime }\mathit{dA}\end{array}$$E32

$$\frac{\partial }{\partial t}A\left[\alpha {\rho }_g\left(h_l^o-\frac{P}{{\rho }_g}\right)+\left(1-\alpha \right){\rho }_l\left(h_g^o-\frac{P}{{\rho }_l}\right)\right]+\frac{\partial }{\partial z}A\left[\alpha {\rho }_g{h}_g^o{v}_g+\left(1-\alpha \right){\rho }_l{h}_l^o{v}_l\right]=q{\prime }_w$$E33

Where

$$h^o=h_k+\frac{1}{2}{v}_k^2+\mathit{gz}$$E34

The subindexes $g$ and $l$ represent the gas and liquid phases, the subindex $i$ represents the property evaluated at the interface, $A$ is the area, $m$ is the average mass transfer per unit of interface area, $P$ is the average pressure, $q$ is heat, $u$ is the velocity at the interface, $w$ is property evaluated at the wall, $\overline{\overline{\sigma }}$ and $\tau$ are wall stresses.

The reader is recommended to refer to reference [17] for further information on this model.

4.3.3 Drift flux model

In earlier two-phase flow models, the fluid was treated as a homogeneous mixture of liquid and vapour, and consequently, only three conservation equations were needed to describe the two-phase flow. In the HEM, the phases are assumed to move at the same velocity, while also experiencing the same temperature. An extension to the mixture model is the drift flux model, in which the relation between the phasic velocities is described through an algebraic equation, thus allowing for a slip between the phases. The Zuber and Findlay [19] drift flux model is probably the most common model that takes into account the relative velocity of the phases. This model uses a drift velocity $v_{\mathit{gj}}$ and a void distribution parameter$C_o$. This model uses the following relationships:

$$j_l=\left(1-\alpha \right)v_l$$E35

$$j_g=\alpha v_g$$E36

$$j=j_l+j_g$$E37

$$v_{\mathit{gj}}=v_g-j=\left(1-\alpha \right)v_g-v_l$$E38

Where $\alpha$ denotes the gas local time fraction, that is, the local void fraction, and $v_k$ is the component, along the axis of the pipe, of the local time-averaged velocity.

We have locally:

$$\alpha v_{\mathit{gj}}=\alpha v_g-\mathit{\alpha j}={\alpha }_j-\mathit{\alpha j}$$E39

Averaging this equation over the total cross-section of the pipe, we obtain the Zuber and Findlay void equation:

$$⟨⟨\alpha ⟩⟩=\frac{⟨⟨j_g⟩⟩}{C_o⟨⟨j⟩⟩+{\stackrel{\sim }{v}}_{\mathit{gj}}}=\frac{J_g}{C_{o}J+{\stackrel{\sim }{v}}_{\mathit{gj}}}$$E40

Where $J_g$ and $J$ are the gas and mixture superficial velocities. In this equation, two quantities appear: (1) The $C_o$ parameter, which accounts for the shape of the $\alpha$ and $\mathrm{j}$ profiles, and (2) the void weighted average of the local drift velocity ${\stackrel{\sim }{v}}_{\mathit{gj}}$, which takes into account the relative velocities between the phases. The symbol $⟨⟨⟩⟩$ defines the average value of a variable $x$ over the cross-sectional area $A$ as

$$⟨⟨x⟩⟩=\frac{1}{A}\underset{A}{\int }x\mathit{dA}$$E41

To calculate $⟨⟨\alpha ⟩⟩$ from this equation is sufficient to know $J_g$, $J$, that is, the volumetric flow rates of each phase and the cross-sectional area, the distribution parameter $C_o$.

$$C_o=\frac{⟨⟨\mathit{\alpha j}⟩⟩}{⟨⟨\alpha ⟩⟩⟨⟨j⟩⟩}$$E42

and the local drift velocity ${\stackrel{\sim }{v}}_{\mathit{gj}}$:

$${\stackrel{\sim }{v}}_{\mathit{gj}}=\frac{⟨⟨\alpha v_{\mathit{gj}}⟩⟩}{⟨⟨\alpha ⟩⟩}$$E43

Zuber and Findlay [19] noticed that $C_o$ and ${\stackrel{\sim }{v}}_{\mathit{gj}}$ are only functions of the flow regime. The authors, therefore, recommended certain values presented in the next table. Nowadays, there are several expressions for ${\stackrel{\sim }{v}}_{\mathit{gj}}$ and $C_o$ available in scientific literature [20].

Note this drift flux model can be applied to homogeneous two-phase flow models to take into account the relative velocities between the gas and liquid phases. However, care must be taken to ensure the distribution parameter $C_o$ and drift velocity ${\stackrel{\sim }{v}}_{\mathit{gj}}$ are applicable to the type of problem and flow regime.

Some models for distribution parameter $C_o$ and drift velocity ${\stackrel{\sim }{v}}_{\mathit{gj}}$ are presented in Table 1.

Flow regime	Distribution parameter $C_o$	Drift velocity
Bubbly	Circular cross-section: $\mathit{\text{Reduced Pressure}}\equiv P_R=\frac{P}{P_c}$ $D>5\mathit{cm}$ $C_o=1-0.5P_R$ $D<5\mathit{cm}$ $P_R<0.5C_o=1.2$ $P_R>0.5C_o=1.4-0.4P_R$	${\stackrel{\sim }{v}}_{\mathit{gj}}=1.41{\left(\frac{\mathit{\sigma g}\mathrm{\Delta }\rho }{{\rho }_l^2}\right)}^{\frac{1}{4}}$
Slug	$C_o=1.2$	${\stackrel{\sim }{v}}_{\mathit{gj}}=0.35{\left(\frac{g\mathrm{\Delta }\mathit{\rho D}}{{\rho }_l}\right)}^{1/4}$
Annular	$C_o=1.0$	${\stackrel{\sim }{v}}_{\mathit{gj}}=23{\left(\frac{{\mu }_l{j}_l}{\mathit{\rho gD}}\right)}^{\frac{1}{2}}\left(\frac{\mathrm{\Delta }\rho }{{\rho }_l}\right)$

Table 1.

Distribution parameter $C_o$ and drift velocity ${\stackrel{\sim }{v}}_{\mathit{gj}}$ [15].

Notice that one of the most important characteristics of multi-phase flow is that it can take on different spatial distributions. This characteristic poses a major challenge to an exact solution treatment [17]. In fact, the majority of the practical two-phase models for conduits, such as the ones presented in this section, were developed for a specific flow regime(s). The following section describes in greater detail the flow regimes for vertical and horizontal flows. Table 2 presents different two-phase flow models based on different assumptions and simplifications [15].

Model designation	Restriction		No. of field equations			No. of interface transfer equations	External constitutive equations
	No	Imposed on	Mass	Momentum	Energy		No	Type
1V1T	3	$v_g=v_l$, $h_l$, $h_G$	1	1	1	0	2	$\tau$,$\overline{q}$
1VS1T		$\frac{v_g}{v_l}$, $h_l$, $h_G$					3	$\tau$,$\overline{q}$, $\frac{v_g}{v_l}=\mathit{\text{slip}}$
1VD1T		$v_g,v_l$, $h_l$, $h_G$					3	$\tau$,$\overline{q,}{v}_r=v_g-v_l$ or $v_{\mathit{gm}}=V_g-v_m$ or $v_{\mathit{gj}}=v_g-v_j$
1VTKTSAT	2	$v_g=v_l$, $h_l$	1	1	2	1	3	$\tau$, $\overline{q}$, $q_g$, $E$
1VSTKTSAT		$\frac{v_g}{v_l}$	1	1	2		4	$\tau$, $\overline{q}$, $q_g$, $\frac{v_g}{v_l}$, $E$
1VDTKTSAT		$v_g-v_l,h_G$	2	1	1		4	$\tau$, $\overline{q}$,$\mathrm{\Gamma }$, $v_r\left(\mathit{\text{or}}{v}_{\mathit{gm}}\mathit{\text{or}}{v}_{\mathit{gj}}\right)$
2V1T		$h_l$, $h_G$	1	2	1		4	${\tau }_l$, ${\tau }_g,\overline{q}$, $M$
1V2T	1	$v_g=v_l$	2	1	2	2	5	${\tau }_l$, $q\prime {\prime }_l$,$q\prime {\prime }_g$, $\mathrm{\Gamma }$, $E$
1VD2T		$v_g-v_l$	2	1	2		6	${\tau }_l$, $q\prime {\prime }_l$,$q\prime {\prime }_g$, $\mathrm{\Gamma }$, $E,v_r\left(\mathit{\text{or}}{v}_{\mathit{gm}}\right)$
2VTKTSAT		$h_l$ or $h_G$	2	2	1		5	${\tau }_l$, ${\tau }_{g,}\overline{q}$, $q_k$, $M$
2VTKTSAT		$h_l$or$h_G$	1	2	2		6	${\tau }_l,{\tau }_g,\overline{q},q_k$, $M$,$E$
2V2T	0	None	2	2	2	3	7	${\tau }_l,{\tau }_g,q_l,q_g,\mathrm{\Gamma }$, $M$,$E$

Table 2.

Practical two-phase flow models [15].

5. Flow regimes

The hypothetical experiment used to explain the boiling curve in the previous subsection describes some of the different regimes a flow can experience during a boiling process. For example, if the coolant undergoes a boiling process, the liquid-gas mixture will present different field structures that give a macroscopic behaviour to the mixture. These structures, also known as regimes or patterns, are classified by visual observation. For that reason, there is no unique set of flow regimes. In thermal-hydraulics analysis, the determination of the flow regime is key, as several variables are dependent on the flow regime, such as the heat transfer coefficient.

In the case of a gas-liquid mixture, gravity plays an important role, as the lighter phase tends to separate from the heavier phase, therefore, the flow regime depends on the orientation of the system. There are different regime classifications for vertical, horizontal and inclined systems.

One classification for vertical and horizontal orientation is presented in Figure 6(a) and (b) respectively.

Figure 6.

Flow regimes for (a) vertical (b) horizontal co-current flows.

The best way to understand each regime is to examine the behaviour of the coolant as it flows through a fuel channel, as shown in Figure 7. The fuel rod is assumed to produce heat, and the flow is assumed to be forced. Normally, one flow pattern transitions into another when the heat flux is changed, or when the mass flux is raised or lowered. For example, subcooled boiling can easily become saturated boiling, and saturated boiling can immediately become bulk boiling. The transition points between these flow regimes are both mass flux and pressure dependent. They depend on the saturation temperature, as well as the coolant channel geometry. However, for a circular coolant channel, these regimes are relatively simple to understand. A graphical depiction of how these flow regimes interact is called a flow regime map.

Figure 7.

Flow regimes for a vertical heated channel.

In some fuel channels of a CANDU-6 reactor, the coolant enters the fuel channel as a subcooled liquid and exits the fuel channel as a two-phase mixture. Under these conditions, the velocities of the liquid and vapour phases are different. Thus, heat transfer coefficients become flow regime dependent.

5.1 Vertical flows

Masterson [7] gives a detailed analysis of flow patterns commonly observed in vertical co-current flows, such as those occurring in Pressurised Water Reactor (PWR) and Boiling Water Reactor (BWR) cores.

Regime I: Subcooled Liquid Flow. When a coolant enters the inlet of a reactor fuel assembly, it is normally a subcooled liquid. The temperature of the coolant stays below the saturation temperature until it progresses further into the core.

Regime II: Bubbly Flow with Subcooled Boiling. Once the temperature of the cladding reaches the saturation temperature, small bubbles begin to form on the surface of the rods. These bubbles detach from the surface and flow into the turbulent core where the bulk fluid temperature is at least several degrees below the saturation temperature ($T_{\mathit{\text{BULK}}}<T_{\mathit{SAT}}$). Any bubbles that form in this way are immediately engulfed by the cooler liquid, and they collapse back into the flow stream. The creation of vapour bubbles at the wall surface is called nucleate boiling and these bubbles do not merge together to form larger bubbles or voids. Hence, the bulk fluid temperature stays below the saturation temperature, and the two-phase mixture tends to be highly turbulent.

Regime III: Bubbly Flow with Saturated Boiling. Bubbly flow with saturated boiling is similar to bubbly flow with subcooled boiling, except that the average temperature of the fluid has now reached the saturation temperature. When this occurs, the bubbles that detach from the surface of the cladding do not immediately collapse when they enter the flow stream. Instead, they either remain intact or combine with other bubbles to form larger bubbles or voids.

Regime IV: Saturated or Bulk Boiling Saturated or bulk boiling is similar to bubbly flow, except that the void fraction now rises from 40 to 60%.

Regime V: Slug Flow. In the slug flow regime, the bubbles that have already formed coalesce into very large bubbles that not only span the entire width of the channel but also remain long and continuous. In reactors, these vapour bubbles can be 4–6 times longer than they are wide. The liquid film on the surface of the rods cannot be agitated anymore, because there is no additional turbulent mixing to increase the heat transfer rate.

Regime VII: Annular Flow. The primary difference between the annular flow regime and the churn flow regime is that the liquid film on the surface of the fuel rods is now moving in the same direction as the vapour. This causes the “churn” in the flow pattern to nearly disappear because the interfacial friction between the two phases is now lower (due to a reduction in the relative velocity difference between the phases). Sometimes this interfacial friction is called interfacial drag.

Regime IX: Pure Vapour Flow. In this final flow regime, the liquid completely disappears and all that remains is pure vapour flow. Hence, this regime is called the pure vapour flow regime. The heat transfer coefficients behave in exactly the same way as they do for single-phase flow. Only in this case, the single phase is the vapour phase.

5.2 Horizontal flows

The flow patterns for co-current horizontal flow are applicable to CANDU reactors. There are several studies that were performed under horizontal and adiabatic (air-water) test sections. However, adiabatic horizontal air-water systems using fuel bundles are limited.

Osamusali et al. [21] presented the generally accepted horizontal flow pattern classifications. These are:

Stratified flow: The stratified flow regime is characterised by the liquid flowing at the bottom of the test section, and the gas phase at the top. This can be further classified: (1) stratified smooth flow. In this case, a smooth gas-liquid interface exists. At high-flow rates, the interface may become wavy. In this case, it is referred to as stratified wavy flow pattern. In fuel bundles such as those used in CANDU reactors, interfacial waves may result from disturbances at the end plates.

Intermittent flow: The intermittent flow regime is characterised by liquid bridging the gap between the gas-liquid interface and the top of the pipe. The liquid bridges are separated by stratified flow zones. The intermittent flow is subdivided into the plug flow regime, occurring at low gas velocities and having a liquid bridge free of gas bubbles, and the slug flow regime, which occurs at higher gas-flow rates and entrains a significant amount of gas bubbles in the liquid bridge. During intermittent flows in rod bundles, the liquid bridges across the elements at the upper part of the channel.

Annular Flow: The annular flow pattern is characterised by the liquid phase flowing around the inner periphery of the pipe and surrounding a core of a fast-flowing gas phase. The gas core may entrain some liquid droplets, and the gas-liquid interface is generally wavy. At low gas-flow rates, the liquid essentially flows as a thick film at the bottom of the pipe with rather unstable waves at the gas-liquid interface, continuously swept up around the pipe periphery, resulting in the wavy-annular flow regime. This eventually leads to the fully developed annular flow regime at higher gas-flow rates, characterised by a continuous liquid film around the inner periphery of the pipe. During annular flows, the rod-bundle elements in the gas core may be covered with very thin liquid films.

Bubbly Flow: The bubbly flow regime is characterised by the void being in the form of discrete bubbles, which are distributed throughout the continuous liquid phase that otherwise fills the pipe section. The bubble concentration is highest at the top of the pipe, especially at lower mass velocities.

Osamusali et al. [21] performed a series of experiments to study the transition of flow patterns in a CANDU-type fuel channel. These authors concluded that the descriptions of two-phase flow patterns occurring in fuel bundles are similar to those observed in pipes. However, Yang [22] observed in a more recent study, that for crept fuel channels the bundle geometry has a strong effect on flow patterns. This effect was pronounced at the transition between stratified and plug/slug flow. For the transition between stratified wavy annular and wavy annular flow, the effect of crept fuel channel was small or non-existent. Bundle misalignment did not show an impact on flow stratification.

5.2.1 Flow regime analysis for CANDU fuel channels

Osamusali et al. [21] developed a generalised flow regime map for a horizontal 37-element bundle, obtained using air-water experimental data at room temperature and near atmospheric pressure. The dimensionless representation is based on dimensionless quantities $X,F,\&T.$ Osamusali et al. [21] presented the following model:

$$X=\sqrt{{\left(\frac{\mathit{dP}}{\mathit{dx}}\right)}_l/{\left(\frac{\mathit{dP}}{\mathit{dx}}\right)}_g}$$E44

The parameter $X$ can be expressed as a function of the superficial velocities:

$$X={\left(\frac{{\mu }_l}{{\mu }_g}\right)}^{0.108}{\left(\frac{{\rho }_l}{{\rho }_g}\right)}^{0.392}{\left(\frac{j_l}{j_g}\right)}^{0.892}$$E45

Where $j_l$ and $j_g$ are the superficial velocities of the gas and liquid phases respectively. The superficial velocity represents the case when the gas of the liquid is assumed to be flowing alone in the channel. For fuel bundles, the superficial velocity is based on the channel cross-sectional area excluding the area occupied by the fuel elements. Using a friction factor developed for a bare 37-element fuel bundle, $f=0.243\mathbf{R}{\mathbf{e}}^{-0.216}$, Eq. (45) can be recast as:

$$X={\left(\frac{{\mu }_l}{{\mu }_g}\right)}^{0.108}{\left(\frac{{\rho }_l}{{\rho }_g}\right)}^{0.5}{\left(\frac{1-x}{x}\right)}^{0.892}$$E46

The dimensionless parameter $\mathbf{F}$ is the modified Froude number, thus given as:

$$\mathbf{F}=\sqrt{\left(\frac{{\rho }_g}{\left({\rho }_l-{\rho }_g\right)g{D}_h}\right)}{j}_g$$E47

$D_h$ is the hydraulic diameter, obtained using the channel cross-sectional area excluding the area occupied by the fuel elements, and the total wetted perimeter of the fuel elements and wall channel. $F$ can be expressed in terms of the total mass flux $G$ and flow quality $x$ as

$$\mathbf{F}=\sqrt{\left(\frac{Gx}{{\rho }_g\left({\rho }_l-{\rho }_g\right)g{D}_h}\right)}$$E48

6. Pressure drop models

Pressure drop is calculated by writing the conservation equations of mass, momentum and energy, and then evaluating the pressure difference. However, in order to solve the basic conservation equations in thermal-hydraulic system codes, additional constitutive equations are required, simply because the number of unknowns is higher than the number of equations. Among them, a constitutive equation is needed to take into account the shear stress.

As mentioned in Section 2, the power generated in the reactor core is removed by the coolant mass flow rate. In a CANDU 6, this mass flow is driven by 4 pumps, which have been dimensioned according to pressure drop calculations and other parameters. Therefore, the design of the reactor plant requires knowledge of the pressure losses across the individual components of the plant. Furthermore, the pressure drop is needed to assess various postulated events and to dimension additional system structures and components needed to support the design requirements of the plant.

In a closed conduit, such as a nuclear fuel channel in a CANDU-6, the pressure drop can be calculated by adding individual pressure drop components as follows:

$$\Delta P_{\mathit{tot}}={\mathit{\Delta P}}_f+\Delta P_{\mathit{acc}}+\Delta P_{\mathit{\text{grav}}}+\Delta P_K$$E49

where $\Delta P_{\mathit{tot}}$ is the total pressure drop, $\Delta P_{\mathit{fr}}$ is the pressure due to the frictional resistance, $\Delta P_{\mathit{acc}}$ is the pressure drop due to acceleration of the flow, $\Delta P_{\mathit{\text{grav}}}$ is the pressure drop due to gravity, and $\Delta P_K$ is the pressure drop due to local flow obstructions.

6.1 Frictional resistance and friction factor

The frictional resistance is the shear stress between the flow and the contacting wall. To account for this shear stress, a non-dimensional friction factor, or the Darcy-Weisbach equation, is commonly used to interrelate the frictional pressure drop to the wall shear stress [23]:

$$f=\frac{\left(\frac{\mathit{dP}}{\mathit{dz}}\right)D}{\frac{1}{2}\rho u^2}$$E50

Where the frictional pressure gradient is negative.

Several correlations and models exist for estimating the friction factor. For example, for turbulent flows in a smooth conduit, the Blasius approximation is widely used [24]:

$$\begin{array}{ccc}f=0.316{\mathbf{Re}}^{-0.25}& \mathit{\text{for}}& 4\times {10}^3<\mathbf{Re}<1\times {10}^5\end{array}$$E51

The Filonenko correlation is also recommended [24]:

$$\begin{array}{ccc}f=\frac{1}{{\left(1.82{log}_{10}\mathbf{R}{\mathbf{e}}_{\mathrm{b}}-1.64\right)}^2}& \mathit{\text{for}}& 4\times {10}^3<\mathbf{Re}<{10}^{12}\end{array}$$E52

However, these correlations only apply to smooth pipes. In reality, most surfaces present some roughness, such as pressure tubes in a CANDU reactor. For that reason, Colebrook combined the smooth wall and fully rough relations in an implicit formula (23):

$$\begin{array}{ccc}f={\left[-2log\left(\frac{ϵ/D}{3.7}\right)+\frac{2.51}{\sqrt{f}\mathbf{Re}}\right]}^{-2}& \mathit{\text{for}}& \mathbf{Re}>4\times {10}^3\end{array}$$E53

However, these models are only valid for isothermal flows or flows without large changes in properties (especially viscosity and density, as these are sensitive to temperature variations).

Furthermore, there are friction factor correlations that were developed for CANDU 37-element fuel bundles [25].

Snoek and Ahmad proposed (as cited in [25]):

$$\begin{array}{ccc}f=0.050\mathbf{R}{\mathbf{e}}^{-0.057}& \mathit{\text{for}}& 108000\le \mathbf{Re}\le 418000\end{array}$$E54

$$Venkat Raj (As cited in [25]) proposed the following correlations for a horizontal 37-element bundle with split-wart spacers, and included the effect of the junctions:

$$f=0.22\mathbf{R}{\mathbf{e}}^{-0.163}\mathit{\text{for}}10000\le \mathbf{Re}\le 140000$$E55

$$f=0.108\mathbf{R}{\mathbf{e}}^{-0.108}\mathit{\text{for}}140000\le \mathbf{Re}\le 50000$$E56

To take into account the non-isothermal (non-iso) nature of these types of flows in calculating the frictional pressure loss, the common approach consists of introducing correction factors into an isothermal friction factor correlation ($f_{\mathit{iso}}$) These correction factors usually take the form of ratios between a fluid bulk property (such as viscosity or density) and the property evaluated at the wall conditions. For example, the following correlations are generally used for supercritical water conditions:

The Kirillov correlation (as cited in [24]):

$$f_{\mathit{\text{non}}-\mathit{iso}}=f_{\mathit{iso}}{\left(\frac{{\rho }_w}{{\rho }_b}\right)}^{0.4}$$E57

Leung and Groeneveld proposed (as cited in [25]):

$$f_{\mathit{\text{non}}-\mathit{iso}}=f_{\mathit{iso}}{\left(\frac{{\mu }_b}{{\mu }_w}\right)}^{-0.28}$$E58

6.2 Acceleration

For a one-dimensional flow in a conduit with axial density variation, the pressure drop due to acceleration can be calculated as [24]:

$$\Delta P_{\mathit{acc}}=G^2\left(\frac{1}{{\rho }_{\mathit{\text{out}}}}-\frac{1}{{\rho }_{\mathit{\text{in}}}}\right)$$E59

6.3 Gravity

This pressure drop component takes into account the gravity force acting on the mass flow. This term, usually called static head, is computed as [24]:

$$\Delta P_{\mathit{\text{grav}}}=\overline{\rho }\mathit{gL}sin\theta$$E60

Where $\overline{\rho }$ is the average density between the inlet and outlet test section, $L$ is the length of the test section, and $\theta$ is the inclination angle.

6.4 Flow obstructions

This component takes into account the obstacles the fluid experiences in a flow, such as endplate spacers in a CANDU fuel channel. This component is generally computed as [24]:

$$\mathrm{\Delta }{P}_K=K\frac{1}{2}\rho u^2=K\frac{G^2}{2\rho }$$E61

Where $K$ is a form loss coefficient usually obtained from experimental data or correlations.

6.5 Two-phase pressure drop

In CANDU-6 reactors the flow is allowed to change phase at the end of some fuel channels. The maximum allowed thermodynamic quality is 4%. Therefore some fuel channels experience two-phase flow conditions. In addition, there are several postulated accidents where two-phase flow conditions are predominant, such as LOCA events. The pressure drop under these two-phase conditions differs greatly from single-phase flow conditions, because the lighter phase tends to travel faster than the heavier phase to satisfy the mass conservation, adding (1) the pressure drop component called acceleration pressure drop$\mathrm{\Delta }{P}_{\mathit{acc}}$, and (2) the pressure losses between the flow and the wall is affected by the presence of a second phase and the flow pattern. In diabatic cases, such as a fuel channel, the heated surface of the element changes the properties near the heated wall, which in turn changes the boundary layer, and thus affecting the friction factor. There are several empirical models to predict the pressure drop under two-phase conditions. Most of them are based on the homogeneous and the separated-flow models (see previous Section 4.3). Most commonly, empirical correlations are used to determine a two-phase friction multiplier. There are several models available, and these vary according to the degree of complexity and might depend on the flow regime.

The two-phase multiplier approach, which is the basis for most of the cited methods, is the generally accepted engineering model to account for the effect of a two-phase mixture in a flow channel. The idea behind this approach is to calculate the pressure drop of one phase (gas or liquid) $\Delta P_L$ first. To determine the two-phase pressure drop $\Delta P_{\mathit{TP}}$, the single-phase pressure drop is multiplied with a two-phase multiplier ${\mathrm{\Phi }}_{2_{l,g}}$ to consider the influence of the second phase. The following four basis methods to define the two-phase multiplier are commonly used: – two methods that assume the liquid or gas phase flowing alone in the flow channel—two methods that assume the entire mixture flowing as liquid or gas only. Applying these basis definitions of the two-phase multipliers, the two-phase friction pressure gradients can be expressed as follows [26]:

$${\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}}_{\mathit{TPF}}=\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}{\mathrm{\Phi }}_l^2$$E62

$${\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}}_{\mathit{TPF}}=\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}{\mathrm{\Phi }}_g^2$$E63

$${\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}}_{\mathit{TPF}}=\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}{\mathrm{\Phi }}_{\mathit{lo}}^2$$E64

$${\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}}_{\mathit{TPF}}=\left\{\frac{\mathrm{\Delta }P}{\mathrm{\Delta }z}\right\}{\mathrm{\Phi }}_{\mathit{go}}^2$$E65

The ${\mathrm{\Phi }}_L^2,{\mathrm{\Phi }}_G^2{\mathrm{\Phi }}_{\mathit{LO}}^2,{\mathrm{\Phi }}_{\mathit{GO}}^2$terms constitute the two-phase pressure drop multiplier to be determined.

6.5.1 Modifications of the single-phase friction factor based on the homogeneous model

In this case, the flow is treated as a pseudo-single-phase fluid, and the friction pressure drop is calculated using a modified friction factor, $f_{\mathit{TP}}$. In this case, the shear stress at the wall is assumed as:

$$\mathrm{\Delta }{P}_{\mathit{TP}}=f_{\mathit{TP}}\frac{L}{D_h}\left(\frac{u}{2{\rho }_m{A}^2}\right)$$E66

Where the $f_{\mathit{TP}}$ is the two-phase friction factor, $\overline{\rho }$ is the mixture density. The friction factor is calculated using friction factor correlations such as Eqs. (53)–(56).

To calculate the $\mathbf{Re}$ of the flow, mixture properties are used. For example, the mixture viscosity can be estimated as:

$$\frac{1}{\mu }=\frac{x}{{\mu }_g}+\frac{1-x}{{\mu }_l}$$E67

or

$$\mu =x{\mu }_g+\left(1-x\right){\mu }_l$$E68

Empirical calculations applying the two-phase multiplier concept based on the separated flow model.

The frictional pressure drop in two-phase flow is generally based on separate flow models, such as the one presented in Section 4.3.2. Recall that in the separated flow model the phases are considered to flow separately in a flow channel, as each phase has its own velocity.

Martinelli and Nelson [27] and Lockhart and Martinelli [28] developed the base for the two-phase friction multiplier approach. Lockhart and Martinelli [28] defined a parameter $\chi$, known as the Martinelli parameter:

$${\chi }^2=\frac{{\mathrm{\Phi }}_g^2}{{\mathrm{\Phi }}_l^2}=\frac{{\left(\mathit{dP}/\mathit{dz}\right)}_l}{{\left(\mathit{dP}/\mathit{dz}\right)}_g}$$E69

The study of two-phase pressure drop has been subject to numerous investigations over the last six decades. Most of the developed methods are based on empirical models.

Finally, the reader is referred to a technical book prepared by an IAEA Coordinated Research Project (CRP), which contains several models for predicting heat transfer coefficients and pressure drop in different flow regimes [25].

Nomenclature

$A$	coefficient $\frac{3.63P}{V{E}_R{\mathrm{\Sigma }}_f}$ (Neutrons/m^2s)
$C_p$	Specific heat capacity (J/kg°C)
$C_o$	distribution parameter (−)
$D$	diameter (m)
$E$	energy transfer between phases (J)
$E_d$	energy deposited locally in the fuel per fission (J)
$E_R$	recoverable energy (J/fission)
$f$	friction factor (—)
$G$	mass flux (kg/m^2s)
$g$	acceleration of gravity (m/s^2)
$h$	specific enthalpy (J/kg)
$\stackrel{\sim }{h}$	heat transfer coefficient (W/m^2°C)
$H_e$	height extrapolated length (m)
$\overline{\overline{I}}$	irreversibility or lost work (J/s)
$\overline{\overline{J}}$	generalised surface source or sink for mass, momentum and energy (—)
$j$	superficial velocity (m/s)
$J_o$	Bessel function of order zero (—)
$K$	form loss coefficient (—)
$k$	thermal conductivity (W/m°C)
$L$	test section length (m)
$M$	momentum transfer between the phases (kg/m s)
$P$	pressure (pa)
${\dot{Q}}_g$	power (W)
$q^{\prime },$	linear heat generation rate (W/m)
$q^{\prime \prime },$	heat flux (W/m^2)
$q^{\prime \prime \prime }$	volumetric heat generation rate (W/m^3)
$\overline{q}$	heat flux between the two-phase and the mixture and the wall (W/m^2)
$r$	radius coordinate (—)
$R_e$	radius extrapolated length (m)
$s_g$	generation of a property per unit of volume and time
$T$	temperature (°C)
$t$	time (s)
$u$	interfacial velocity (m/s)
$v$	velocity (m/s)
$V$	volume (m^3)
$v_{\mathit{gj}}$	drift velocity (m/s)
$X$	dimensionless parameter (—)
$x$	quality (—)
$z$	axial coordinate (—)

Greek letters

$\alpha$	void fraction (—)
$\mathrm{\Gamma }$	net vapour volumetric production rate resulting from phase change (kg/m^3 s)
$\gamma$	angle (degrees)
$\mathit{\delta w}$	phase evapouration rate (kg/s)
$ϵ$	roughness (m)
$\mathrm{\Theta }$	neutron flux (Neutrons/m^2s)
$\theta$	angle (rad)
${\chi }^2$	two-phase pressure drop parameter “Martinelli” (—)
${\mathrm{\Phi }}^2$	two-phase pressure drop multiplier (—)
$\mu$	viscosity (kg/m s)
$\phi$	property per unit of volume (1/m^3)
$\rho$	density (kg/m^3)
${\mathrm{\Sigma }}_f$	macroscopic fission cross-section (cm^−1)
$\overline{\overline{\sigma }}$	normal stress component (N/m^2)
$\tau$	shear stress component (N/m^2)

Non-dimensional numbers

$\mathbf{F}$	modified Froude number $\mathbf{F}\equiv \sqrt{\left(\frac{{\rho }_g}{\left({\rho }_l-{\rho }_g\right)g{D}_h}\right)}{j}_g$ (—)
$\mathbf{Nu}$	Nusselt number, $\equiv \frac{\mathit{hD}}{k}$ (—)
$\mathbf{Pr}$	Prandtl number, $\equiv \frac{C_p\mu }{k}$ (—)
$\mathbf{Re}$	Reynolds number, $\equiv \frac{G{D}_e}{\mu }$ (—)
${\mathit{P}}_{\mathit{r}}$	reduced pressure, $\equiv \frac{P}{P_c}$ (—)

Mathematical symbols

$⟨⟨⟩⟩$	Cross section average (—)
$\nabla$	Del (also known as Nabla) operator (—)
$\overrightarrow{}$	vector (—)

Subscripts or superscripts

acc	acceleration
b	bulk
c	critical
$\mathit{\text{centre}}$	property evaluated at the centre
e	extrapolated length
f	friction factor
g	gas phase
grav	gravity
in	inlet
iso	isothermal
K	form loss
k	phase ($g$=gas, $l$=liquid)
l	liquid phase
non-iso	non-isothermal
out	outlet
R	reduced
tot	total
TPF	two-phase flow
w	wall

Acronyms and abbreviations widely used in text and list of references

CANDU	CANada Deuterium Uranium
CRP	Coordinated Research Project
D2O	deuterium
DNB	Departure Nucleate Boiling
EFP	End-Flux-Peaking
HEM	Homogeneous Equilibrium Model
IAEA	International Atomic Energy Agency
LOCA	loss-of-coolant accident
NPP	Nuclear Power Plant