two-phase flow transport and heat transfer

5

One-/two-phase flow transport and heat transfer 5.1 Tools required for flow transport modelling 5.1.a

Introduction

The modelling of one- and two-phase flow transport is a very challenging task for many different reasons. First, even if a flow is at steady-state conditions, some small fluctuations around the mean value of the flow properties are always present. For one-phase flows, these fluctuations originate from turbulence. For two-phase flows, these fluctuations also come from the strongly time-dependent character of the relative spatial distribution of the two phases. Furthermore, since nuclear reactors are large and complex systems, it is still impossible with the available computing power of today to resolve all the meaningful scales for the encountered flow situations. As a result, the description of one- and two-phase flow dynamics in nuclear reactors is carried out at a macroscopic level, both with respect to space and time. In concrete terms, this means that the local balance equations, that were introduced in Section 2.3 and in Section 2.4 of Chapter 2, for expressing the conservation of mass and momentum and of energy/enthalpy, respectively, have to be time-averaged and space-averaged on proper time bins and elementary volumes, respectively, so that the complexity of the problem at hand is reduced. This is also equivalent to saying that the high-frequency and small-scales phenomena are filtered out, thus simplifying the problem. As will be seen in Section 5.2 of this chapter, this time- and space-averaging procedure introduces new unknowns to the problem being modelled, resulting in a larger number of unknowns than the number of available equations. These unknowns actually correspond to the effect of the high-frequency and small-scale phenomena onto the macroscopic balance equations. In order to close the system of equations and in addition to the equation of state of the fluid, other relationships are needed. Such relationships either describe the transfer of mass, momentum and energy/enthalpy at the interfaces between the phases (jump conditions), or the transfer of momentum and energy/enthalpy at the interface between any of the two phases and the walls the phases might be in contact with. These relationships strongly depend on the flow configuration or morphology, usually referred to as flow regimes. Finally, a so-called topological law is required to compensate for the fact that the information about the relative distribution of the phases was lost in the averaging process. Modelling of Nuclear Reactor Multi-physics. https://doi.org/10.1016/B978-0-12-815069-6.00005-2 Copyright © 2020 Elsevier Inc. All rights reserved.

251

252

Modelling of Nuclear Reactor Multi-physics

In this section, the different flow regimes for two-phase vertical flows usually encountered in light-water reactors (LWRs) are highlighted. Thereafter, the mathematical tools necessary to carry out the time- and space-averaging of the local conservation equations are derived.

5.1.b

Two-phase flow regimes

As mentioned above, the macroscopic models used for expressing the conservation of mass, momentum and energy of two-phase flows rely on a number of constitutive relationships, usually determined experimentally, as well as on the equation of state of the fluid. These additional relationships are needed to close the system of equations, together with the boundary and initial conditions. The constitutive relationships depend on the flow configuration or morphology, usually referred to as flow regimes or flow patterns. In vertical flows, the main flow regimes, schematically represented in Fig. 5.1, are (Todreas and Kazimi, 1993)
Single-phase liquid. Only the liquid phase is present.
Bubbly flow. This flow regime is characterized by the presence of dispersed vapour bubbles in a continuous liquid phase. The bubbles can be of variable size and shape.
Slug flow. This flow regime is characterized by the presence of gas plugs (or large bubbles) separated by liquid slugs. The liquid film that surrounds the gas plugs usually moves downward. In addition, several small bubbles may be dispersed within the liquid.
Churn flow. This flow regime has essentially the same features as for slug flow. The main difference is a much more pronounced ‘chaotic’ behaviour.
Annular flow. This flow regime is characterized by the presence of a continuous core of gas surrounded by an annulus of liquid. In case of a sufficiently high gas flow, the core of gas may carry liquid droplets, and the flow is then called annulardispersed flow.
Drop flow. There is not any longer any annulus of liquid. The core of gas still carries liquid droplets.
Single-phase vapour. Only the vapour phase is present. Depending on the relative ‘properties’ of the liquid phase as compared to the vapour phase, a given flow regime exists. The mapping of the various flow regimes as a function of the relative ‘properties’ of the liquid phase as compared to the vapour phase is referred to as flow regime map. By ‘properties’, some measure of the fluxes of each of the two phases is considered, although other parameters and variables might be included in a flow regime map. Such flow regime maps can only be determined experimentally. Many flow regime maps exist, depending on the chosen ‘properties’, the facility in which the experiments are performed, the configuration of the experiments, etc. Various flow

Chapter 5
One-/two-phase flow transport and heat transfer 253

FIGURE 5.1 Main flow regimes in vertical flows. Derived from Todreas, N.E., Kazimi, M.S., 1993. Nuclear Systems I: Thermal Hydraulic Fundamentals. Taylor & Francis, Levittown, USA.

regime maps can be found, and the interested reader is referred to the existing literature on this subject, such as Hewitt and Roberts (1969), Taitel and Dukler (1977), Taitel et al. (1980). It should be emphasized that a flow regime map is specific to the conditions and experimental set-up the map was derived for. The use of a flow regime map for

254

Modelling of Nuclear Reactor Multi-physics

modelling a nuclear reactor fuel assembly is only valid if the conditions of the experiments are indeed present in the considered nuclear fuel assembly. Finally, one should observe that the ‘properties’ of the liquid and vapour phases, on which a flow regime map is established, are not known in advance: the goal of the thermal-hydraulic calculations is to determine the properties of the fluid, by solving macroscopic balance equations. Such equations involve experimental correlations for some of the terms. These correlations are strongly dependent on the flow regime, which itself depends on the ‘properties’ of the liquid and vapour phases. An iterative solution procedure is thus required.

5.1.c

Mathematical tools

Tools necessary to properly integrate the local balance equations with respect to time and space are introduced in order to derive time- and space-averaged macroscopic balance equations. Although this is the most common procedure used in nuclear thermal-hydraulics, it should be mentioned that other techniques using, e.g., ensemble averaging can also be utilized to obtain balance equations. Ensemble averaging relies on the basic idea that having access to many realizations or events, the corresponding expected value is the average over all those events (Berry, 2003). Balance equations for those averages can be correspondingly derived. Time- and space-averaging techniques can be seen as a subset of ensemble averaging (Drew and Passman, 1999). Ensemble averaging techniques are nevertheless superior to time- and space-averaging techniques, since ensemble averaging does not rely on the necessity to have control volumes containing a sufficient quantity of a given component of the fluid to give meaningful ensemble averages, as a classical space-averaging technique would otherwise require. In the following, only time- and space-averaging techniques are presented. The interested reader is referred to the literature existing on ensemble-averaging techniques.

Descriptors in two-phase flow dynamics For the purpose of time- and space-averaging, different descriptors are introduced. It should be pointed out that these averaged properties are still time- and spacedependent, since the time bins do not encompass the whole duration of the possible reactor transient and the volume elements do not encompass the whole system. In order to better describe the two phases, a phase density function is introduced and defined as:
At time t; Xk ðr; tÞ ¼

1 if r ˛ phase k 0 if r ; phase k

(5.1)

If at position r and time t, the phase is of type k (k ¼ l for the liquid phase, or k ¼ v for the vapour or gaseous phase), then the phase density function is equal to one. Otherwise, the phase density function is equal to zero. Different types of averages can be defined: space-averages, time-averages and spaceand time-averages (Todreas and Kazimi, 1993). The averaging with respect to space can

Chapter 5
One-/two-phase flow transport and heat transfer 255

be performed either on a specific phase k or on both phases. Concerning the space averages, the averaging can be performed on typically a cross-sectional surface area A of the flow or on a volume element V of the flow. In the following, both will be generically denoted as a domain D, where D thus represents either A or V. The instantaneous (two-phase) space-average of a quantity c(r, t) on a domain D is defined as: hci ¼

Z

1 D

cðr; tÞd n r

(5.2)

D

where


n ¼ 2 if DhA n ¼ 3 if DhV

(5.3)

The instantaneous (one-phase) space-average of the quantity c(r, t) on the domain Dk occupied by the phase k is defined as: hck i ¼

1 Dk

Z cðr; tÞd n r

(5.4)

Dk

Since on Dk, one has c(r, t) ¼ ck(r, t), Eq. (5.4) can also be written as: hck ik ¼

1 Dk

Z

ck ðr; tÞd n r

(5.5)

Dk

Both notations are equivalent, i.e. hck ik ¼ hck i

(5.6)

The local time-average of a quantity c(r, t) on a time interval Dt is defined as: c¼

1 Dt

Z

tþDt 2 tDt 2

cðr; tÞdt

(5.7)

If Dtk represents the cumulative time during which the phase k is present at a given point r in the time interval Dt, another local time-average can be defined as: 1 c¼ Dtk

Z

tþ

Dtk 2

Dtk 2

cðr; tÞdt

(5.8)

t

As opposed to the time-average estimated according to Eq. (5.7), the time-average given by Eq. (5.8) is estimated on the cumulative time Dtk , i.e. on the time during which the phase k is present at a given point r in the total time interval Dt. After introducing the various space- and time-averaging procedures above, some properties associated to those averaging procedures are detailed hereafter.

Time-averaging properties Following the previous definitions, with respect to the phase density function Xk(r, t), one has

256

Modelling of Nuclear Reactor Multi-physics


For the local time-average of Xk(r, t) on Dt, usually denoted as ak : Z

1 Dt

a k ¼ Xk ¼

tþDt 2

Xk ðr; tÞdt

tDt 2

(5.9)


For the local time-average of Xk(r, t) on Dtk , usually denoted as ak : ak ¼ Xk ¼

Z

1 Dtk

Dtk 2

tþ

Xk ðr; tÞdt

Dtk 2

(5.10)

t

Since Xk(r, t) ¼ 0 during times when the phase k does not occupy the position r, it can be seen that: ak ¼

1 Dtk

Z

Dtk 2

tþ

Xk ðr; tÞdt ¼ 1

Dtk 2

(5.11)

t

and ak ¼

1 Dt

Z

tþDt 2 tDt 2

Xk ðr; tÞdt ¼

1 Dt

Z

tþ

Dtk 2

Dtk 2

Xk ðr; tÞdt ¼

t

Dtk Dt

(5.12)

Space-averaging properties Following the previous definitions, the instantaneous (two-phase) space-average of Xk(r, t), usually denoted as hak i, is given as: hak i ¼ hXk i ¼

1 D

Z

Xk ðr; tÞd n r ¼

Dk D

(5.13)

D

The space-average of the gaseous phase density function, usually denoted hai, is thus defined as: hai ¼ hav i ¼

Dv Dv ¼ D Dl þ Dv

(5.14)

and corresponds to the void fraction. The volume-average of the liquid density function is then given as: hal i ¼ 1  hai

(5.15)

Correspondingly, with respect to ck(r, t), one has
For the instantaneous (two-phase) space-average of ck(r, t): hck i ¼

1 D

Z

ck ðr; tÞd n r ¼ Dk

1 D

Z

ck ðr; tÞXk ðr; tÞd n r

(5.16)

D


For the instantaneous (one-phase) space-average of ck(r, t): hck ik ¼

1 Dk

Z

ck ðr; tÞd n r ¼ Dk

1 Dk

Z

ck ðr; tÞXk ðr; tÞd n r D

(5.17)

Chapter 5
One-/two-phase flow transport and heat transfer 257

from which one notices, using Eqs. (5.13) and (5.16), that: hck i ¼ hck ik

Dk ¼ hck ikhak i D

(5.18)

Space-time-averaging properties Double-averaging is considered hereafter. The space-time-average of Xk(r, t), usually denoted as hak i , is defined as: 1 hak i ¼ hXk i ¼ Dt

Z

0

tþDt 2

@1 D

tDt 2

Z

1

Xk ðr; tÞd rAdt n

The space-time average of ck(r, t) is given by: 1 hck i ¼ Dt

Z

tþDt 2

tDt 2

(5.19)

D

0

@1 D

1

Z

ck ðr; tÞd rAdt n

(5.20)

D

Because of Eq. (5.18), this can be rewritten as: hck i ¼ hck ikhak i

(5.21)

If the instantaneous values of hck ik and hak i can be split into a time-averaged component and a fluctuating part as: hck ik ¼ hck ik þ hck i0k with hck i0k ¼ 0

(5.22)

0 0 hak i ¼ hak i þ hak i with hak i ¼ 0

(5.23)

hck i ¼ hck ik hak i þ hck i0k hak i0

(5.24)

then Eq. (5.21) becomes:

The space-time average of ck is thus usually different from the product between the (one-phase) space-time average of ck and the space-time average of ak . The instantaneous space averages can be time-averaged, and the local time averages can be space-averaged. We will also assume that the space- and time-averaging operators are commutative (Todreas and Kazimi, 1993), i.e. hci ¼

1 Dt

Z

tþDt 2 tDt 2

0 @1 D

Z

1 cðr; tÞd n rAdt ¼

D

1 D

Z D

0 B1 @ Dt

Z

1 tþDt 2

tDt 2

C cðr; tÞdt Ad n r ¼ hci

(5.25)

and one thus has:   hak i ¼ ak

(5.26)

258

Modelling of Nuclear Reactor Multi-physics

Using the property of commutativity, one could also verify that: hak ihck ik 1 ¼ Dt

1 ¼ D

Z

Z D

¼

1 D

Z

0 tþDt 2 tDt 2

BDk 1 @ D Dk

0 B1 @ Dt

Z

Z

1 1 C ck ðr; tÞd rAdt ¼ Dt

Z

n

Dk

1 tþDt 2 tDt 2

1 C ck ðr; tÞXk ðr; tÞdt Ad n r ¼ D

tþDt 2 tDt 2

Z D

0 @1 D

1

Z

ck ðr; tÞXk ðr; tÞd rAdt n

D

0 BDtk 1 @ Dt Dtk

Z

Dtk 2

tþ

Dtk 2

t

1

(5.27)

C ck ðr; tÞdt Ad n r

D E ak ck d n r ¼ ak ck

D

5.2 Derivation of the space- and time-averaged conservation equations for flow transport 5.2.a

Introduction

In the following, the local balance equations derived in Section 2.3 and in Section 2.4 of Chapter 2, for expressing the conservation of mass and momentum and of energy/ enthalpy, respectively, are first space-averaged and then time-averaged. The equations derived in this section would be identical to the ones obtained by first time-averaging the local balance equations and then space-averaging them. The local balance conservation equations can be written in the following generic form: vðrf Þ ðr; tÞ þ V $ ðrf vÞðr; tÞ ¼ ðrqÞðr; tÞ  V $ Fðr; tÞ vt

(5.28)

vðrfÞ ðr; tÞ þ V $ ðrf5vÞðr; tÞ ¼ ðrqÞðr; tÞ  V $ Fðr; tÞ vt

(5.29)

if f is a scalar, and

if f is a vector, where the different quantities appearing in the above equations are given in Table 5.1. It has to be noted that the energy conservation equation can be replaced by the enthalpy conservation, following the derivation already introduced in Section 2.4.b of Chapter 2. The two formulations are completely equivalent.

5.2.b

Space-averaging of the local conservation equations

Using a control volume V that encompasses both the fluid and solid structures (typically the nuclear fuel rods or spacer grids) as represented in Fig. 5.2, the local conservation equations given by Eqs. (5.28) and (5.29) can be volume-averaged on the volume Vk

Chapter 5
One-/two-phase flow transport and heat transfer 259

Table 5.1 Definition of the variables appearing on the right-hand side of Eqs. (5.28) and (5.29) for the conservations of mass, momentum and energy. Conservation of mass

Conservation of momentum

Conservation of energy

f or f F

1 0

v (sPI)

e q00 (sPI)$v

q or q

0

g

q000 þ g:v r

FIGURE 5.2 Schematic representation of a control volume encompassing solids (in dark grey) and two phases (with the liquid phase represented in light grey). The interface between the two phases is indicated by the dashed lines.

occupied by the phase k within the control volume V and then divided by V. This reads as: 1 V

Z Vk

vðrf Þ 3 1 d rþ vt V

Z V $ ðrf vÞd 3 r ¼

Z

1 V

Vk

rqd 3 r 

1 V

Vk

Z V $ Fd 3 r

(5.30)

Vk

if f is a scalar, and 1 V

Z Vk

vðrfÞ 3 1 d rþ vt V

Z V $ ðrf5vÞd 3 r ¼ Vk

1 V

Z rqd 3 r  Vk

1 V

Z V $ Fd 3 r

(5.31)

Vk

if f is a vector. In the above two equations, the explicit space- and time-dependence was dropped for the sake of clarity of the equations. Nevertheless, it has to be emphasized that such equations are still space-dependent, since the control volume chosen for the integration is only a part of the entire computational domain, i.e. the above balance equations need to be written for each of these control volumes.

260

Modelling of Nuclear Reactor Multi-physics

Using the general transport theorem (see Eq. (2.61) of Chapter 2), as well as the generalized Gauss theorem (i.e. Eq. (1.32) of Chapter 1), Eq. (5.30) can be rearranged into: 2 v 61 4 vt V 1 ¼ V

Z

3 7 1 rf d 3 r5  V

Vk

Z

2

I

Z

61 rf v S $ Nd 2 r þ V $ 4 V Z

61 rqd 3 r  V $ 4 V

Vk

7 1 rf vd 3 r5 þ V

Vk

Sk

2

3

3

Z

7 1 Fd 3 r5  V

Vk

Z rf v $ Nd 2 r þ Ski

Z rf v $ Nd 2 r Skw

(5.32)

Z

1 F $ Nd 2 r  V

Ski

1 V

F $ Nd 2 r Skw

where Sk represents the surface of the volume Vk occupied by the phase k within the control volume V. This outer surface can actually be decomposed into the surface being common with the outer surface of the control volume V, the interface Skw between the phase k and the solid walls and the interface Ski between the two phases. Since the control volume is a fixed control volume and since the solid walls are also fixed, on both of which one has vS ¼ 0, one concludes that: 1 V

I

rf v S $ Nd 2 r ¼

1 V

Sk

Z

rf v S $ Nd 2 r

(5.33)

Ski

Furthermore, since the flow is necessarily tangential to the solid walls, one has v$N ¼ 0 on the solid walls. Consequently, Eq. (5.32) simplifies into: 2 v 61 4 vt V

¼

1 V

Z

3

2

7 61 rf d 3 r5 þ V $ 4 V

Vk

Z

Z

3 7 rf vd 3 r5

Vk

2 61 rqd 3 r  V $ 4 V

Vk

Z

(5.34)

3 7 1 Fd 3 r5  V

Vk

Z ½rf ðv  v S Þ þ F $ Nd 2 r 

1 V

Ski

Z F $ Nd 2 r Skw

Likewise, using the general transport theorem (see Eq. (2.62) of Chapter 2), as well as the generalized Gauss theorem (i.e. Eq. (1.33) of Chapter 1), Eq. (5.31) can be rearranged into: 2 v 61 4 vt V þ

1 V

Z

7 1 rfd 3 r5  V

I

Vk

2 61 ðrf5v S Þ $ Nd 2 r þ V $ 4 V

ðrf5vÞT $ Nd 2 r þ Z Vk

2 61 rqd 3 r  V $ 4 V

Z

3 7 rf5vd 3 r5

Vk

Sk

Z Ski

1 ¼ V

3

1 V

Z ðrf5vÞ $ Nd 2 r Skw

Z Vk

3

7 1 Fd 3 r5  V

Z FT $ Nd 2 r  Ski

(5.35)

1 V

Z FT $ Nd 2 r Skw

where Sk represents the surface of the volume Vk occupied by the phase k within the control volume V. This outer surface can actually be decomposed into the surface being

Chapter 5
One-/two-phase flow transport and heat transfer 261

common with the outer surface of the control volume V, the interface Skw between the phase k and the solid walls and the interface Ski between the two phases. Since the control volume is a fixed control volume and since the solid walls are also fixed, on both of which one has vS ¼ 0, one concludes that: 1 V

I

ðrf5v S Þ $ Nd 2 r ¼

1 V

Sk

Z

ðrf5v S Þ $ Nd 2 r

(5.36)

Ski

Furthermore, one has: ðrf5vÞT $ N ¼ ðv5rfÞ $ N ¼ vðrf $ NÞ

(5.37)

For the case of the momentum conservation equation for which f ¼ v, one has v$N ¼ 0 on the solid walls since the flow is necessarily tangential to the walls. This thus means that on the solid walls, one has, based on Eq. (5.37), ðrf5vÞT $ N ¼ 0. Consequently, Eq. (5.35) simplifies into: 2 v 61 4 vt V Z

1 þ V

Z

3

2

7 61 rfd 3 r5 þ V $ 4 V

Vk

Z

3 7 1 rf5vd 3 r5 ¼ V

Vk

Z

2 61 rqd 3 r  V $ 4 V

Vk

  1 ðrf5vÞT þ rf5v S  FT $ Nd 2 r  V

Ski

Z

3 7 Fd 3 r5

Vk

(5.38)

Z F $ Nd r T

2

Skw

For the case of the momentum conservation equation for which f ¼ v and F ¼ (sPI), one has ðrf5vÞT ¼ rf5v and FT ¼ F (since the stress tensor s is symmetric, as indicated by Eq. (2.77) of Chapter 2). One thus finally obtains: 2 v 61 4 vt V 1  V

Z

Z

3

2

7 61 rfd 3 r5 þ V $ 4 V

Vk

Z

3 7 1 rf5vd 3 r5 ¼ V

Vk

1 ½rf5ðv  v S Þ þ F $ Nd r  V

2 61 rqd 3 r  V $ 4 V

Vk

Z

Z Vk

3 7 Fd 3 r5 (5.39)

F $ Nd r

2

Ski

Z

2

Skw

Since the control volume V might include solids, one defines a parameter, called the porosity, giving the ratio between the volume occupied by the fluid and the control volume. Since the fluid might contain both liquid and vapour, the porosity is thus defined as: b¼

Vl þ Vv V

(5.40)

As such, the described procedure for deriving the space- and time-averaged conservation equations for mass, momentum and energy is usually referred to as the porous media approach. Using the porosity, one notices that for any quantity c one can write: 1 V

Z cd3 r ¼ Vk

Vl þ Vv Vk 1 V Vl þ Vv Vk

Z cd3 r ¼ bhak ihck ik Vk

(5.41)

262

Modelling of Nuclear Reactor Multi-physics

where the last equality was obtained using Eqs. (5.13) and (5.17). Based on the relationship given by Eq. (5.41), the volume-averaged conservation equations can be finally written as: v ½bhak ihrk fk ik þ V $ ½bhak ihrk fk v k ik ¼ bhak ihrk qk ik  V $ ½bhak ihFk ik vt Z Z 1 1 2 ½rk fk ðv k  v S Þ þ Fk  $ Nd r  Fk $ Nd 2 r  V V Ski

(5.42)

Skw

if f is a scalar, and v ½bhak ihrk f k ik þ V $ ½bhak ihrk f k 5v k ik ¼ bhak ihrk qk ik  V $ ½bhak ihFk ik vt Z Z 1 1 ½rk f k 5ðv k  v S Þ þ Fk  $ Nd 2 r  Fk $ Nd 2 r  V V Ski

(5.43)

Skw

if f is a vector.

5.2.c

Time-averaging of the space-averaged conservation equations

The volume-averaged equations given by Eqs. (5.42) and (5.43) can be time-averaged on a time interval Dt. Following the notations introduced in Section 5.1.c, this reads as: v ½bhak ihrk fk ik þ V $ ½bhak ihrk fk v k ik ¼ bhak ihrk qk ik  V $ ½bhak ihFk ik vt Z Z 1 1 ½rk fk ðv k  v S Þ þ Fk  $ Nd 2 r  Fk $ Nd 2 r  V V Ski

(5.44)

Skw

if f is a scalar, and v ½bhak ihrk f k ik þ V $ ½bhak ihrk f k 5v k ik ¼ bhak ihrk qk ik  V $ ½bhak ihFk ik vt Z Z 1 1  ½rk f k 5ðv k  v S Þ  Fk  $ Nd 2 r  Fk $ Nd 2 r V V Ski

(5.45)

Skw

if f is a vector. Noticing that: V$c ¼V$c

(5.46)

for c being either a vector or a tensor, and 2 vc v 61 ðtÞ ¼ 4 vt vt Dt

Z

    Dt Dt Z tþDt2 c tþ c t 1 vc 0 0 vc 2 2 7 cðt 0 Þdt 0 5 ¼ ¼ ðt Þdt ¼ ðtÞ Dt Dt tDt2 vt vt 3

tþDt 2

tDt 2

(5.47)

Chapter 5
One-/two-phase flow transport and heat transfer 263

for c being a scalar, or 2 vc v 61 ðtÞ ¼ 4 vt vt Dt ¼

Z

    Dt Dt Z tþDt2 c t c tþ 1 vc 0 0 2 2 7 ¼ ðt Þdt cðt 0 Þdt 0 5 ¼ Dt Dt Dt t 2 vt 3

tþDt 2

tDt 2

(5.48)

vc ðtÞ vt

for c being a vector, Eqs. (5.44) and (5.45) can be transformed, using Eq. (5.27), into: 2 3 3 2 2 * * * * +3 + + + v6 6 7 7 4b ak ðrk fk Þ 5 þ V $ 4b ak ðrk fk v k Þ 5 ¼ b ak ðrk qk Þ  V $ 4b ak Fk 5 vt Z

1  V

Ski

(5.49)

Z

1 ½rk fk ðv k  v S Þ þ Fk  $ Nd 2 r  V

Fk $ Nd 2 r Skw

if f is a scalar, and 2

3 2 * +3 + + * v6 7 7 6 4b ak ðrk f k Þ 5 þ V $ 4b ak ðrk f k 5v k Þ 5 ¼ b ak ðrk qk Þ  V $ 4b ak Fk 5 vt 1  V

*

Z

+

3

2

*

1 ½rk f k 5ðv k  v S Þ þ Fk  $ Nd r  V

Ski

(5.50)

Z Fk $ Nd r

2

2

Skw

if f is a vector. Eqs. (5.49) and (5.50) represent the space- and time-averaged conservation equations, in which the temporal variations on time scales smaller than the time interval Dt and the spatial variations on scales smaller than V1/3 were filtered out. Due to the time- and space-averaging, new quantities have appeared in the conservation equations, quantities associated with transfers at the interfaces and at the solid walls.

5.2.d

Equations to be solved

Using the expressions for the different variables given by Table 5.1, Eq. (5.49) represents the generic conservation equation to be used for mass and energy and Eq. (5.50) represents the generic conservation equation to be used for momentum. Each of these equations can be written for each of the two phases, and one thus obtains six equations. Before turning to the necessary approximations required to solve this set of equations, some further transformations are necessary so that one obtains a set of equations containing the same type of unknowns in all equations. One first ‘de-couples’ the two-phase averages of products between the time-averaged phase density function and any time-averaged phasic quantity (either in scalar form ck or in vector/tensor form ck) into products between the two-phase average of the

264

Modelling of Nuclear Reactor Multi-physics

_

_

time-averaged phase density function and a corresponding averaged quantity ck or ck , respectively, as (Obry, 1996): *

+ ak ck

¼ ak



+

* _ ck

¼

or ak ck



_ ak ck _

(5.51) _

As can be seen from this expression, the average quantities ck and ck are defined as: *

_

ck

+

*

1R ak ck ak ck d 3 r VV _ or ck ¼ ¼ ¼ 1R 3 ak d r ak VV

+

ak ck

ak

1R ak ck d 3 r VV ¼ 1R ak d 3 r VV

(5.52)

and can thus be seen as a two-phase volume average of the time-averaged phasic quantity ck or ck using the local time-averaged phase density function as a weighting function. Likewise, one also needs to ‘de-couple’ the two-phase averages of products between the time-averaged phase density function, any time-averaged phasic quantity (in scalar form ck), and any time-averaged phasic quantity (in scalar form dk or in vector/tensor form dk). The two-phase averages of those products are replaced by products between the two-phase average of the time-averaged phase density function, the average quan_ c tities r as defined above, and a corresponding averaged quantity dc k or dk , respeck

tively, as (Obry, 1996): +

* ak ðrk dk Þ

* +

_ c ck ¼ ak rk dk or ak ðrk dk Þ ¼ ak r_k d

(5.53)

c As can be seen from this expression, the average quantities dc k and dk are defined as: +

*

ak ðrk dk Þ dc k ¼



1R 1R ak ðrk dk Þ d 3 r ak ðrk dk Þ d 3 r VV VV ¼ ¼ _ 1R rk R ak rk d 3 r 3 ak d r VV V V

_ ak rk

or *

+ ak ðrk dk Þ

ck ¼ d



where Eq. (5.52) was used.

_ ak rk

1R 1R ak ðrk dk Þ d 3 r ak ðrk dk Þ d 3 r VV VV ¼ ¼ _ 1R rk R ak rk d 3 r ak d 3 r VV V V

(5.54)

Chapter 5
One-/two-phase flow transport and heat transfer 265

Based on Eqs. (5.51) and (5.53), the space- and time-averaged conservation equations become:
For the conservation of mass: 2 3 2 3



Z v4 1 _ _ b ak rk 5 þ V $ 4b ak rk c vk 5 ¼ ½  rk ðv k  v S Þ $ Nd 2 r vt V

(5.55)

Ski


For the conservation of momentum: 2 2 2 3 2 3 3 3









_ v4 _ _ _ _ b ak rk c v k 5 þ V $ 4b ak rk ðv kd 5v k Þ 5 ¼ b ak rk g þ V $ 4b ak sk 5  V4b ak Pk 5 vt þ

1 V

Z ½  rk v k 5ðv k  v S Þ þ sk  Pk I $ Nd 2 r þ

1 V

Ski

Z ½sk  Pk I $ Nd 2 r Skw

(5.56)

where the following identity was used: 2 2 3 3 3







_ _ _ _ 4 4 4 5 5 V $ b ak sk  b ak Pk I ¼ V $ b ak sk  V $ b ak Pk I5 2

2

2 2 2 3 3 3 3









_ _ _ _ _ ¼ V $ 4b ak sk 5  b ak Pk V $ I  V4b ak Pk 5 $ I ¼ V $ 4b ak sk 5  V4b ak Pk 5 (5.57)


For the conservation of energy: 2 2 3 3







_ v4 _ _ _ 4 5 d b ak rk ebk þ V $ b ak rk ðek v k Þ 5 ¼ b ak q000 þ b ak rk c vk $ g k vt 3 2 3 3 2





_ h h 6 00 7 V $ 4b ak qk 5 þ V $ 4b ak ðsk $v k Þ 5  V $ 4b ak ðPk v k Þ 5 2

1 þ V

Z Ski



 rk ek ðv k  v S Þ þ ðsk  Pk IÞ $ v k 

q000 k



1 $ Nd r þ V

Z

2

Skw

  ðsk  Pk IÞ $ v k  q00k $ Nd 2 r

(5.58)

266

Modelling of Nuclear Reactor Multi-physics

which, since vk$N ¼ 0 on the solid walls, can be rearranged into:

3 2 2 3







_ v4 _ _ _ 000 5 b ak rk ebk 5 þ V $ 4b ak rk ðed vk $ g ¼ b ak qk þ b ak rk c k vk Þ vt

3 2 3 3 2





_ h h 7 6 V $ 4b ak q00k 5 þ V $ 4b ak ðsk $v k Þ 5  V $ 4b ak ðPk v k Þ 5 2

þ

1 V

Z



 1 2  rk ek ðv k  v S Þ þ ðsk  Pk IÞ $ v k  q000 k $ Nd r þ V

Ski

Z

(5.59)

  sk $ v k  q00k $ Nd 2 r

Skw

As can be seen in the above equations, the following cross-terms have appeared: h h 5v k Þ , ðed ðv kd k v k Þ , ðsk $v k Þ and ðPk v k Þ because of the averaging introduced by Eqs. (5.51) and (5.53). These terms can be treated in the following manner. We first split each of the following local phasic quantities into a time average, as defined by Eq. (5.10), and a deviation from such an average (Ishii and Ibiki, 2011): rk ¼ rk þ r0k

(5.60)

v k ¼ v k þ v 0k

(5.61)

ek ¼ ek þ ek0

(5.62)

sk ¼ sk þ s0k

(5.63)

Pk ¼ Pk þ Pk0

(5.64)

Taking the time average, as defined by Eq. (5.10), on Eqs. (5.60)e(5.64) result in: r0k ¼ 0

(5.65)

v 0k ¼ 0

(5.66)

ek0 ¼ 0

(5.67)

s0k ¼ 0

(5.68)

Pk0 ¼ 0

(5.69)

Based on Eqs. (5.61)e(5.64), one has:  v k 5v k ¼ v k þ

v 0k

   0 5 v k þ v k ¼ v k 5v k þ v 0k 5v 0k þ v k 5v 0k þ v 0k 5v k

(5.70)

Chapter 5
One-/two-phase flow transport and heat transfer 267

   0 0 ek v k ¼ ek þ ek0 v k þ v 0k ¼ ek v k þ ek0 v 0k þ ek v k þ ek v k 1 0

0 sk $ v k ¼ @sk þ

1

s0k A $ @v k

0

(5.71)

þ

v 0k A

10

¼ sk $ v k þ s0k $ v k þ sk $ v 0k þ s0k $ v 0k

(5.72)

1

0 0 0 0 Pk v k ¼ @Pk þ Pk0 A@v k þ v 0k A ¼ Pk v k þ Pk v k þ Pk v k þ Pk v k

(5.73)

Because of Eqs. (5.66)e(5.69), Eqs. (5.70)e(5.73) reduce to: v k 5v k ¼ v k 5v k þ v 0k 5v 0k

(5.74)

ek v k ¼ ek v k þ ek0 v 0k

(5.75)

sk $ v k ¼ sk $ v k þ s0k $ v 0k

(5.76)

Pk v k ¼ Pk v k þ Pk0 v 0k

(5.77)

Usually, the following terms are neglected (Analytis, 2003): s0k $ v 0k z0

(5.78)

Pk0 v 0k z0

(5.79)

and

so that one obtains for Eqs. (5.76) and (5.77): sk $ v k z sk $ v k

(5.80)

Pk v k zPk v k

(5.81)

In addition, the fluctuations in density are negligible, i.e. r0k z0, so that Eq. (5.60) simply becomes (Analytis, 2003): rk zrk

(5.82)

Using Eqs. (5.74), (5.75), (5.80) and (5.81) together with Eq. (5.82), the expressions for h h the cross-terms ðv kd 5v k Þ , ðed k v k Þ , ðsk $v k Þ and ðPk v k Þ can be rewritten as: +

*

*

ak rk ðv k 5v k Þ ðv kd 5v k Þ ¼



_ ak rk

+

*

ak rk ðv k 5v k Þ z



_ ak rk



ak rk v k 5v k ¼



_ ak rk



*

+

ak rk þ

v 0k 5v 0k



_ ak rk



+

(5.83)

268

Modelling of Nuclear Reactor Multi-physics *

+

* ak rk ðek v k Þ ðed k vk Þ ¼



+

*

ak rk ðek v k Þ z

_ ak rk



h ðsk $v k Þ ¼

ak rk ek v k ¼

_ ak rk

*

+

*

h ðPk v k Þ ¼

z

_ ak rk

(5.84)

+





*

ak ðPk v k Þ





+

(5.85)

ak +





ak sk $ v k



ak

 ak rk ek0 v 0k

þ



_ ak rk

ak ðsk $ v k Þ

*

*

+





+

ak Pk v k z





ak

(5.86)

ak

Defining the spatial correlation coefficients Cvk ;vk , Cek ;vk , Csk ;vk and CPk ;vk such that: *



ak rk v k 5v k



+

¼ Cvk ;vk ak

+

* ak rk ek v k *



ak sk $ v k *









ak Pk v k

0

_ rk

1

@ vc A c k 5v k



_ c ¼ Cek ;vk ak rk ec k v k

+

+





(5.87)

(5.88)

_

¼ Csk ;vk ak sk $ vc k

(5.89)



_ ¼ CPk ;vk ak Pk vc k

(5.90)

Eqs. (5.83)e(5.86) become: * 0

1

 ak rk v 0k 5v 0k

Aþ c ðv kd 5v k Þ z Cvk ;vk @ vc k 5v k

* c c ðed k v k Þ z Cek ;v k e k v k þ



_ ak rk

 ak rk ek0 v 0k

_ ak rk

+

(5.91)

+

(5.92)

Chapter 5
One-/two-phase flow transport and heat transfer 269

h _ ðsk $ v k Þ z Csk ;vk sk $ vc k

(5.93)

h _ ðPk v k Þ zCPk ;vk Pk vc k

(5.94)

In macroscopic modelling tools as the ones used for nuclear reactors, one finally assumes that all spatial correlation coefficients are equal to unity (Delhaye, 2008): Cvk ;vk z1

(5.95)

Cek ;vk z1

(5.96)

Csk ;vk z1

(5.97)

CPk ;vk z1

(5.98)

The above assumptions are sometimes referred to as topological laws. Topological laws represent assumptions on the phasic spatial correlation coefficients and on the phasic pressures that are necessary to close the system of equations we intend to derive. The purpose of the topological laws is to artificially introduce in the space-/timeaveraged balance equations some information about the local spatial distributions that were lost during the spatial averaging. Equating the spatial correlation coefficients according to Eqs. (5.95)e(5.98) represents the simplest (but necessary) possible assumption on the phasic spatial correlations. The topological laws with respect to pressure will also be detailed in Section 5.3.a. Using those assumptions in Eqs. (5.91) and (5.94) finally leads to: * 0

1

 ak rk v 0k 5v 0k

Aþ c ðv kd 5v k Þ z @ vc k 5v k

* c c ðed k vk Þ z e k v k þ



_ ak rk

 ak rk ek0 v 0k

_ ak rk

+

(5.99)

+

(5.100)

h _ ðsk $ v k Þ z sk $ vc k

(5.101)

h _ ðPk v k Þ zPk vc k

(5.102)

Before giving the final form of the conservation equations, the momentum conservation equation can be further simplified as described hereafter. By introducing

270

Modelling of Nuclear Reactor Multi-physics

the pressure fluctuations p0k as the difference between the instantaneous value of the _

pressure Pk and its average value Pk : _

p0k ¼ Pk  Pk

(5.103)

one obtains: 1 V

Z ½  rk v k 5ðv k  v S Þ þ sk  Pk I $ Nd 2 r þ Ski

¼

1 V

1 V

Z ½sk  Pk I $ Nd 2 r Skw

Z



 1  rk v k 5ðv k  v S Þ þ sk  p0k I $ Nd 2 r  V

Ski

þ

1 V

Z

1 z V

_

Pk I $ Nd 2 r Ski



 1 sk  p0k I $ Nd 2 r  V

Skw

Z

Z

Z

_

Pk I $ Nd 2 r (5.104)

Skw



 1  rk v k 5ðv k  v S Þ þ sk  p0k I $ Nd 2 r þ V

Ski

Z

  sk  p0k I $ Nd 2 r

Skw

2

3

Z 61 Z 7 _ 1 6 7 Nd 2 r þ Nd 2 r 7 Pk I $ 6 4V 5 V Ski

Skw

In virtue of the generalized Gauss theorem (i.e. Eq. (1.32) of Chapter 1) applied for c(r) ¼ 1, the last term on the right-hand side of the previous equation can be rewritten as (Analytis, 2003): 2

3

61 _ 6 Pk I $ 6 4V

Z Nd 2 r þ Ski

1 V

Z Skw



Z 7 _ _ _ 1 7 Nd 2 r 7 ¼ Pk I $ V d 3 r ¼ Pk bVhak i ¼ Pk bV ak 5 V

(5.105)

Vk

where the last equality was obtained noticing that Vc ¼ Vc

(5.106)

and using Eq. (5.25). One thus gets: 1 V

Z ½  rk v k 5ðv k  v S Þ þ sk  Pk I $ Nd 2 r þ Ski

1 z V

1 V

Z ½sk  Pk I $ Nd 2 r Skw

Z Ski



 rk v k 5ðv k  v S Þ þ sk 



p0k I

1 $ Nd r þ V

Z

2



  _ sk  p0k I $ Nd 2 r þ Pk bV ak

(5.107)

Skw

Consequently, the final form of the conservation equations is given as follows:

Chapter 5
One-/two-phase flow transport and heat transfer 271


For mass, by: 3 2 3 2



v4 _ 5 _ v k 5 ¼ Mki b ak rk þ V $ 4b ak rk c vt

(5.108)


For momentum, by: 2 2 3 0 13



v4 _ _ A5 c b ak rk c v k 5 þ V $ 4b ak rk @ vc k 5v k vt 0 2 13





_ _ _ t A5 @ 4  b ak VPk þ Fki þ Fkw sk þ sk zb ak rk g þ V $ b ak

(5.109)


For energy, by: 2 2 3 3







_ v4 _ _ _ 5 c b ak rk ebk 5 þ V $ 4b ak rk ec vk $ g q000 þ b ak rk c k v k zb ak k vt 0 13 2 3 3 2





_ _ _ B 00 6 t 00 C7 5  V $ 4b ak Pk vc 5 V $ 4b ak @qk þ qk A5 þ V $ 4b ak sk $ vc k k 2

(5.110)

þEki þ Ekw

where the different quantities appearing in the above equations are defined as: * ak rk stk ¼ 

v 0k 5v 0k



+

z



ak *

qk ¼



ek0 v 0k

+



Mki ¼ 1 V

Z

+ _

¼ rk





ak rk ek0 v 0k z





*





ak

Fki ¼

 ak rk v 0k 5v 0k

0d 0  v k 5v k

(5.111)

ak

ak rk t 00

*



+ _

¼ rk

d  ek0 v 0k

(5.112)

ak 1 V

Z ½  rk ðv k  v S Þ $ Nd 2 r

(5.113)

  rk v k 5ðv k  v S Þ þ sk  p0k I $ Nd 2 r

(5.114)

Ski

Ski

Fkw ¼

1 V

Z Skw



 sk  p0k I $ Nd 2 r

(5.115)

272

Modelling of Nuclear Reactor Multi-physics

Eki ¼

1 V

Z



  rk ek ðv k  v S Þ þ ðsk  Pk IÞ $ v k  q00k $ Nd 2 r

(5.116)

Ski

Ekw ¼

Z

1 V

  sk $ v k  q00k $ Nd 2 r

(5.117)

Skw 00

where stk represents a stress tensor due to fluctuations/turbulence, qtk represents the heat flux due to fluctuations/turbulence; Mki is the mass transfer through the liquid/ vapour interface due to phase change; Fki represents the transfer of momentum through the liquid/vapour interface due to phase change and the forces due to stresses and pressure fluctuations; Eki represents the transfer of energy through the liquid/vapour interface due to phase change, the work of the forces due to pressure/ stresses and temperature gradients; Fkw represents the forces at the solid walls due to stresses and pressure fluctuations; and Ekw represents the work of the forces at the solid walls due to stresses and the transfer of energy at the solid walls. At the liquid/vapour interface, the following so-called jump interface conditions are fulfilled: Mli þ Mvi ¼ 0

(5.118)

Fli þ Fvi ¼ 0

(5.119)

Eli þ Evi ¼ 0

(5.120)

Eqs. (5.118), (5.119) and (5.120), respectively, express the fact that the mass, momentum and energy transferred from the liquid phase to the vapour phase are equal to the ones received by the vapour phase from the liquid phase. This assumes that there is no source/sink of mass, momentum and energy, respectively, at the interfaces. Such an assumption is valid for mass, whereas for momentum and energy, this assumption is equivalent to neglecting surface tension compared to the other forces in presence (Delhaye, 2008). In the conservation equations given by Eqs. (5.108)e(5.110), the averaged quantities are defined according to Eqs. (5.51) and (5.53), i.e. *

+

1R ak rk d 3 r VV

¼ ¼ 1R ak d 3 r ak VV ak rk

_ rk

(5.121)

+

*

1R ak sk d 3 r VV ¼ ¼ 1R ak d 3 r ak VV ak sk

_

sk

(5.122)

Chapter 5
One-/two-phase flow transport and heat transfer 273 +

* ak Pk _ Pk

¼

ak

1R ak Pk d 3 r VV ¼ 1R ak d 3 r VV

(5.123)

1R ak ðrk v k Þ d 3 r VV ¼ 1R ak rk d 3 r VV

(5.124)

1R ak ðrk ek Þ d 3 r VV ¼ 1R ak rk d 3 r VV

(5.125)

and +

* ak ðrk v k Þ vc k ¼



_ ak rk +

* ak ðrk ek Þ ec k ¼



_ ak rk

As earlier mentioned, the conservation equation for energy is completely equivalent to a conservation equation written for enthalpy. The conservation equation of enthalpy could be derived in the same manner as for energy from the local balance equation for enthalpy. Nevertheless, such a procedure is much more tedious. Furthermore, some difficulties arise from the fact that the local enthalpy is defined as h ¼ u þ P=r with _

e ¼ u þ v2/2, and thus the definition of the average enthalpy hk following Eq. (5.53) _

_

leads to some incompatibilities with the average density rk , the average pressure Pk , vk , i.e. and the average velocity c h _ v k Þ2 Pk ðc ðv k Þ2 Pk _ þ þ _ h k ¼ ek  s ek  2 rk 2 rk _

(5.126)

Instead, one defines an average enthalpy as (Ishii and Hibiki, 2011): .

_ h k ¼ ek



v k Þ2 ðc 2

_

þ

Pk _

rk

(5.127) _

For obtaining a balance equation for enthalpy, Eq. (5.109) is first multiplied by v k in a scalar manner as: 8 2 8 2 39 0 1 39






0 0 139 1



= _ _ v  b ak @VPk A $ c v k þ Fki $ c v k þ Fkw $ c vk þ V $ 4b ak @sk þ stk A5 $ c : ; k 2

(5.128)

274

Modelling of Nuclear Reactor Multi-physics

Noticing that: 8 8 2 2 39 0 139



<


< = = 1



_ vk b ak rk c 2

12

@c vk A v $ c vk þ vt

2

2 3

v4 _ 5 b ak rk vt

8 0 2 2 3 139



= 1< _ _ v k 4V $ c vk 5 þ c v k $ 4V5@b ak rk c v k A5 $ c v þ b ak rk c ; k 2: 8 0 2 2 13 39



= 1< _ _ 5 $c c v k $ 4V $ @b ak rk c v k A5 þ b ak rk c v k $ 4V5 vc v þ k ; k 2: 0 8 2 39




_ vk b ak rk c 2

12

@c vk A v $ c vk þ vt

2

2 3

v4 _ b ak rk 5 vt

8 0 2 2 3 139



= 1< _ _ v k 4V $ c vk 5 þ c v k $ 4V $ @b ak rk c v k A5 $ c v þ b ak rk c ; k 2: 8 0 2 2 13 39



= 1< _ _ 5 $c c v k $ 4V5@b ak rk c v k A5 þ b ak rk c v k $ 4V5 vc v þ k ; k 2: 2

0

12 3

6 @c vk A 6

v6 _ 6 ¼ 6b ak rk vt 6 2 6 4

0

12

7 @c 2 3 vk A 7

7 v 7 4b ak r_k 5 7þ 7 vt 2 7 5

9 8 2 3 3 2



= 1< v _ _ v k 4V $ c vk 5  c v k $ 4b ak rk 5 þ c v k $ Mki $ c v þ b ak rk c ; k 2: vt 8 2 0 139

= 1< 4 _ v k $ @b ak rk c v k A5 $ c v þ V c ; k 2:

(5.129)

Chapter 5
One-/two-phase flow transport and heat transfer 275

where the mass conservation equation (i.e. Eq. 5.108) was used in the last equality, as well as the following identity: 0

2 2 0 13 3 13





_ _ _ 5 ¼ V4 c c v k $ 4V5@b ak rk c v k A5 þ b ak rk c v k $ 4V5 vc v k $ @b ak rk c v k A5 k 2

(5.130)

one finds that: 8 8 2 39 0 139 2



= = <
0

12 3

6 @c vk A 6

6 v6 _ ¼ 6b ak rk vt 6 2 6 4 2

0

7 8 2 0 2 3 139 0 12 7

= 7 1<

1 _ _ 7 v k 4V $ c v k 5 þ V4 c v k $ @b ak rk c v k A5 $ c v þ Mki @ c vk A b ak rk c 7þ ; k 7 2: 2 7 5

12 3

6 @c vk A 6

6 v6 _ ¼ 6b ak rk vt 6 2 6 4

2

0

12

3

7 7 6 @c 0 12 vk A 7 7 6

7 7 1 6 _ 7 7 6 c v k 7 þ Mki @ c vk A 7 þ V $ 6b ak rk 7 7 2 6 2 7 7 6 5 5 4 (5.131)

Eq. (5.128) can thus be rewritten as: 2

0

12 3

6 @c vk A 6

v6 _ 6 6b ak rk vt 6 2 6 4

2

0

12

3

7 7 6 @c 0 12 vk A 7 7 6

7 7 1 6 _ 7 7 6 c v k 7 þ Mki @ c vk A 7 þ V $ 6b ak rk 7 7 2 6 2 7 7 6 5 5 4 (5.132)



_ vk þ ¼ b ak rk g $ c

8 < :

2

0



_ 4 V $ b ak @sk

0 139 1

= _ v  b ak @VPk A $ c vk þ stk A5 $ c ; k

v k þ Fkw $ c vk þFki $ c

or, using the identity: 0 139

= _ v V $ 4b ak @sk þ stk A5 $ c ; k : 9 2 8 0 0 1 13

= <

_ _ tA t A5 @ @ 4 : V5 c v k  b ak vk ¼ V $ b ak sk þ sk $ c sk þ sk ; : 8 <

2

(5.133)

276

Modelling of Nuclear Reactor Multi-physics

as: 2

0

12 3

6 @c vk A 6

6 v6 _ 6b ak rk vt 6 2 6 4

2

0

12

3

7 7 6 @c vk A 7 7 6

7 7 6 _ 7 7 6 c vk 7 7 þ V $ 6b ak rk 7 7 6 2 7 7 6 5 5 4

9 2 8 0 0 1 13



= <

_ _ _ tA t A5 @ @ 4 : V5 c v k þ V $ b ak v k  b ak vk sk þ sk $ c sk þ sk ¼ b ak rk g $ c ; :

(5.134)

0 1 0 12

_ 1 v k  Mki @ c v k A þ Fki $ c v k þ Fkw $ c vk b ak @VPk A $ c 2

This equation expresses the conservation of kinetic energy. Subtracting this equation to the conservation equation for the energy (i.e. Eq. 5.110), and using the definition of the average enthalpy as introduced in Eq. (5.127), one obtains: 2 3 2 3







. _ _ v4 _ . _ 5 b ak rk hk  b ak Pk 5 þ V $ 4b ak rk hk vc Pk vc k  b ak k vt 13 0 2 0 3 2 13 2







_ _ 00 _ C 7 B 6 5 þ 4b ak @sk þ st A5: V5 c zb ak qk000  V $ 4b ak @q00k þ qtk A5  V $ 4b ak stk $ vc vk k k 0 3 1 0 12



_ _ 1 5 þ b ak @VPk A $ c v k þ Mki @ c v k A  Fki $ c v k  Fkw $ c v k þ Eki þ Ekw V $ 4b ak Pk vc k 2 2

(5.135)

which simplifies into: 2 2 3 3



. . v4 _ _ 5 b ak rk hk 5 þ V $ 4b ak rk hk vc k vt 2 13 0 2 0 3 2 13







_ _ 00 6 C7 B 00 5 þ 4b ak @s_k þ st A5: V5 c zb ak q000 vk @qk þ qtk A5  V $ 4b ak stk $ vc k k  V $ 4b a k k 0 3 1 0 12 2



_ _ v4 1 v k þ Mki @ c v k A  Fki $ c v k  Fkw $ c v k þ Eki þ Ekw þ b ak Pk 5 þ b ak @VPk A $ c vt 2 (5.136)

Chapter 5
One-/two-phase flow transport and heat transfer 277

5.3 Flow models In the following, different flow models are presented. Those are the two-fluid model based on the set of equations derived in the above section; a mixture model in which the two phases are in thermal equilibrium but move at different velocities; and a homogeneous equilibrium model in which both phases are in thermodynamic equilibrium and move at the same velocity. Those models represent the main categories of flow models. Although many variations of such models exist, the following section will entirely focus on those models. The interested reader is referred to the existing literature on the other models. In the averaging process described above, the local balance equations are replaced by phasic balance equations, i.e. the local balance equations for mass, momentum and energy/enthalpy were transformed in equations volume-averaged on each of the phases (and time-averaged). This results in losing information about the relative distribution of the phases between each other and can be seen as the phases coexisting with their respective phase-averaged properties at all spatial points through each of the control volumes. This approach is thus sometimes referred to as the interpenetrating media formulation.

5.3.a

Two-fluid model

In the two-fluid model, the mass, momentum and energy/enthalpy conservation equations are time- and space-averaged on each of the two phases. The information about the spatial variations of the phases within the control volumes as well as about the temporal variations of the phases within the time bins is lost, i.e. small scales and highfrequency phenomena are filtered out. The equations that one needs to solve are the time- and space-averaged conservation equations earlier derived and given:
For mass, by Eq. (5.108), i.e. 2 2 3 3



v4 _ _ b ak rk 5 þ V $ 4b ak rk c v k 5 ¼ Mki vt

(5.137)


For momentum, by Eq. (5.109), i.e. 2 3 0 13 2



v4 _ _ A5 c v k 5 þ V $ 4b ak rk @ vc b ak rk c k 5v k vt 0 2 13





_ _ _ zb ak rk g þ V $ 4b ak @sk þ stk A5  b ak VPk þ Fki þ Fkw

(5.138)

278

Modelling of Nuclear Reactor Multi-physics


For energy, by Eq. (5.110), i.e. 2 2 3 3







_ v4 _ _ _ 5 c b ak rk ebk 5 þ V $ 4b ak rk ec vk $ g q000 rk c k v k zb ak k þ b ak vt 13 0 2 3 3 2





_ _ 00 _ C7 B 6 5  V $ 4b ak Pk vc 5 V $ 4b ak @q00k þ qtk A5 þ V $ 4b ak sk $ vc k k 2

(5.139)

þEki þ Ekw

or alternatively for enthalpy, by Eq. (5.136), i.e. 2 2 3 3



. v4 _ . 5 _ 5 b ak rk hk þ V $ 4b ak rk hk vc k vt 13 0 2 0 3 2 13 2







_ _ _ B 00 6 000 t 00 C7 t t 5 þ 4b ak @sk þ s A5: V5 c zb ak qk  V $ 4b ak @qk þ qk A5  V $ 4b ak sk $ vc vk k k 0 3 1 0 12 2



_ _ v4 1 b ak Pk 5 þ b ak @VPk A $ c v k þ Mki @ c v k A  Fki $ c v k  Fkw $ c v k þ Eki þ Ekw þ vt 2 (5.140)

These equations are valid on each of the two phases. It thus means that one has six equations. The model above is thus also referred to as a six-equation model. The unknowns to be determined are _


The average phasic density rk (as defined by Eq. 5.121).
The average phasic velocity c v k (as defined by Eq. 5.124).
The average phasic (total) energy ebk (as defined by Eq. 5.125) or alternatively the .

average phasic enthalpy hk (as defined by Eq. 5.127). These unknowns, sometimes referred to as primary variables, must be determined on each of the two phases, and thus amount to six quantities. When looking at the conservation equations, one notices that these equations contain additional unknown quantities, called secondary variables. It should be emphasized that the choice of the primary and secondary variables is not unique. When implementing a fluid solver, the primary variables will be solved explicitly for, whereas the secondary variables will only be determined implicitly, based on additional relationships that need to be introduced. In the present book, the secondary variables are as follows:
Terms related to transfers at the liquid/vapour interfaces: the mass transfer Mki through the liquid/vapour interface due to phase change (as defined by Eq. 5.113),

Chapter 5
One-/two-phase flow transport and heat transfer 279

the transfer of momentum Fki through the liquid/vapour interface due to phase change and the forces due to stresses and pressure fluctuations (as defined by Eq. 5.114), the transfer of energy Eki through the liquid/vapour interface due to phase change, the work of the forces due to pressure/stresses and temperature gradients (as defined by Eq. 5.116).
Terms related to transfers between each of the two phases, respectively, and the solid walls: the forces Fkw at the solid walls due to stresses and pressure fluctuations (as defined by Eq. 5.115) and the work Ekw of the forces at the solid walls due to stresses and the transfer of energy at the solid walls (as defined by Eq. 5.117). These unknowns must be determined on each of the two phases, and thus amount to 10 quantities. In addition, additional terms are present in the space- and time-averaged conservation equations. Those are as follows: D E
The average phase density function ak . _
The average pressure Pk . _


The average stress tensor sk .
The stress tensor stk due to fluctuations/turbulence. _


The average heat flux q00k . 00


The heat flux qtk due to fluctuations/turbulence. _


The average volumetric heat source q000 k . These unknowns must be determined on each of the two phases, and thus amount to 14 quantities. In summary, one has six equations and 30 unknown quantities. In order to have a problem that can be solved, 24 additional equations/relations are necessary in order to get as many unknowns as equations/relations, as explained below (Obry, 1996; Delhaye, 2008). The jump interface conditions, given by Eqs. (5.118)e(5.120) and recalled hereafter: Mli þ Mvi ¼ 0

(5.141)

Fli þ Fvi ¼ 0

(5.142)

Eli þ Evi ¼ 0

(5.143)

give three additional relationships. In virtue of the definition of the average phase density function (see Eq. 5.13), one also has:

al

þ

av

¼1

(5.144)

280

Modelling of Nuclear Reactor Multi-physics

This gives one additional relationship. The equation of state of the fluid gives a relationship between three state variables on each of the two phases, as: rk hrk ðPk ; hk Þ

(5.145)

Although the equation of state relates the local instantaneous variables rk , Pk and hk, it is generally assumed that the equation of state also holds for the space- and time_

_

.

averaged variables rk , Pk and hk , i.e. one can also write (Obry, 1996): 0

1

_ _ _ rk hrk @Pk

.

; hk A

(5.146)

This equation of state for the averaged state variables can be written on each of the two phases, thus giving two additional relationships. In addition to the equation of state, other constitutive equations relating some of the fluid properties are available. A constitutive equation gives a relationship between the heat flux and other quantities, as: q00k hq00k ðPk ; hk ; Tk Þ

(5.147)

with Tk being the phasic temperature, linked to the state variables via the equation of state that can also be written as: Tk hTk ðPk ; hk Þ

(5.148)

Such a constitutive equation can for instance be based on Fourier’s law of conduction, already introduced in Chapter 2, and which is recalled hereafter: q00k ¼  kk ðTk ÞVTk

(5.149)

Furthermore, a constitutive equation gives a relationship between the stress tensor and other quantities, respectively, as: sk hsk ðPk ; hk ; v k Þ

(5.150)

Such a constitutive equation can, for instance, be based on the assumption of a Newtonian fluid, according to which one has for the stress tensor:   2 sk ¼ mk V5v k þ ðV5v k ÞT  mk IðV $ v k Þ 3

(5.151)

with mk being the dynamic viscosity of the phase k (known from the fluid characteristics), and where the superscript T represents the transpose. As before, the above constitutive relationships hold for local instantaneous variables, and it is assumed that the same constitutive relationships can also be applied to time- and space-averaged variables, i.e. 0 _ q00k

_ hq00k

_ @Pk

1 .

fk A ; hk ; T

(5.152)

Chapter 5
One-/two-phase flow transport and heat transfer 281

and 0 _ sk

_ hsk

_ @Pk

with

1 .

vk A ; hk ; c

0

(5.153)

1

_ fk hTk @Pk T

.

; hk A

(5.154)

These constitutive equations can be written on each of the two phases, thus giving four additional relationships. The two phases are not necessarily in thermal equilibrium (or more exactly in thermodynamic equilibrium). This means that in thermal non-equilibrium, one has: _

_

Tl sTv

(5.155)

The interfacial correlations used would need in such a case to account for such a _

thermal non-equilibrium and an extra correlation giving a relationship between Tl and _

Tv would also be required. In the case of thermal equilibrium, one has on the other hand: _

_

Tl ¼ Tv ¼ Tsat

(5.156)

The two phases have thus the same temperature, equal to the saturation temperature Tsat. In thermal equilibrium, the equation of state for the fluid only involves relationships between two thermodynamic variables. _ Usually, the volumetric heat source q000 is a given input to the problem, and this k provides two additional inputs/relationships. A modelling of the turbulent terms (i.e. the turbulent stress tensor stk and the tur00 bulent heat flux qtk ) is required for each of the two phases. This thus means that the modelling of four terms is required. The turbulent stress tensor stk is usually given, by analogy with the stress tensor (i.e. Eq. 5.151), by   stk ¼ mtk V5v k þ ðV5v k ÞT

(5.157)

mtk

where is a parameter having the meaning of a turbulent dynamic viscosity. Such a parameter cannot be known from the fluid characteristics, and a model is thus required 00 for it. Likewise, the turbulent heat flux qtk is usually given, by analogy with Fourier’s law of conduction (i.e. Eq. 5.149), by: 00

qtk ¼  ck Vhk

(5.158)

where ck is a parameter having the meaning of a turbulent ‘diffusivity’ for the enthalpy. Such a parameter cannot be known from the fluid characteristics, and a model is thus required for it.

282

Modelling of Nuclear Reactor Multi-physics

A modelling of the transfer of mass/momentum/energy at the liquid/vapour interface (i.e. Mki , Fki and Eki ) and of the transfer of momentum/energy at the solid walls (i.e. Fkw and Ekw ) is required to close the system of equations. Because of Eqs. (5.141)e(5.143), the modelling of the transfer at the liquid/vapour interfaces is only required on one of the two phases. For the transfer at the solid walls, modelling is required for each of the two phases. This thus means that the modelling of seven terms is required. Many different models exist for each of these terms, based on first-principle (mechanistic), semi-empirical or completely empirical approaches. Presenting the various models is beyond the scope of this book. The interested reader is instead referred to the extensive literature existing on this subject. Summarizing the number of unknowns and the number of available relationships, one notices that one has in total 30 unknowns and 29 available equations/relationships if thermal equilibrium is assumed. Therefore, an additional relationship is required to close the system of equations. The ‘missing’ relationships come from the fact that in the averaging process, some information was lost about the relative distribution of the two phases and their interface structure and dynamics. As a consequence, an additional topological law is required to close the system of equations. Typically, this additional topological law relates the pressure distributions. Such a law depends very much on the type of flows being studied and their configuration. Basically, two types of topological laws with respect to pressure are used. They both relate the pressures in each phase and the pressure at the interface. The first type of law is to assume pressure equilibrium, i.e. that in the direction perpendicular to the main flow transport, all pressures are equal: _

_

Pl ¼ Pv ¼ Pi

(5.159)

where Pi represents the pressure at the liquid/vapour interface. In such a case, the pressure fluctuations p0k appearing in the interfacial terms (i.e. in Eq. 5.114) and in the wall terms (i.e. in Eq. 5.115) are identically equal to zero. This approximation, which is most often used, is usually acceptable for vertical flows, although, as will be explained in Section 5.4, such an approximation leads to a non-hyperbolic character of the equations to be solved, possibly giving rise to numerical instabilities. The other possible topological law with respect to pressure is to have pressure non-equilibrium, i.e. that in the direction perpendicular to the main flow transport, the pressures are not equal: _

_

Pl s Pv sPi

(5.160)

The difficulty in letting the pressures be different is that a relationship between the phasic pressures and the interfacial pressure is thus necessary. There are unfortunately only very few situations where the phasic and interfacial pressures can be determined. For instance, in case of horizontal stratified flows having a simple interfacial geometry, such a relationship can be easily derived assuming a hydrostatic pressure distribution in the direction perpendicular to the flow direction (Delhaye, 2008). Some attempts have also been made to include surface tension in the momentum balance equation at the

Chapter 5
One-/two-phase flow transport and heat transfer 283

liquid/vapour interfaces. The difficulty in such approaches is to specify the surface tension in complicated non-steady-state interfacial geometries. In most two-fluid models, the phasic pressures and interfacial pressures are taken equal in the direction perpendicular to the direction of the main flow transport.

5.3.b

Mixture models with specified drift velocities

One possible simplification of the two-fluid model is to sum the phasic equations in order to obtain a system of equations expressing the conservation of mass, momentum and energy/enthalpy for the mixture. The main advantage of this approach lies with the fact that the interfacial transfer terms disappear due to the jump interface conditions, as will be demonstrated hereafter. The modelling of those terms, rather complex and difficult in the two-fluid model, is not required in mixture models, thus greatly simplifying the implementation of such models. Nevertheless, it will be demonstrated that such a methodology introduces a drift velocity between the phases as a new unknown in the momentum conservation equation, for which a ‘closure’ relationship is required. Summing up the mass conservation equations (i.e. Eq. 5.137) of the two-fluid model and using the jump interface condition (i.e. Eq. 5.141) lead to: 2 2 3 3







v4 _ _ 5 _ _ b al rl þ b av rv þ V $ 4b al rl vbl þ b av rv c vv 5 ¼ 0 vt

(5.161)

This equation can be written as: b

v r þ V $ ðbGm Þ ¼ 0 vt m

(5.162)

where the mixture density is given as:



_ _ _ _ rm ¼ al rl þ av rv ¼ ð1  haiÞrl þ hairv

(5.163)

and the mixture mass flux is given as:



_ _ _ _ v v ¼ ð1  haiÞrl vbl þ hairv c vv Gm ¼ al rl vbl þ av rv c

(5.164)

In the above two equations, Eq. (5.144) was used and the average density function for   the vapour phase av is simply denoted as hai. Summing up the momentum conservation equations (i.e. Eq. 5.138) of the two-fluid model and using the jump interface condition (i.e. Eq. 5.142) lead to, using Eqs. (5.163) and (5.164): 2 0 13 X _ v 4b ak rk @ vc A5 c b Gm þ V $ k 5v k vt k¼l;v

¼ brm g þ V $

X k¼l;v

0 13 2 3

X _ X 4b ak @s_k þ st A5  b 4 ak VPk 5 þ Fkw k 2

k¼l;v

k¼l;v

(5.165)

284

Modelling of Nuclear Reactor Multi-physics

Assuming uniform pressure distribution between the phases and using Eq. (5.144), Eq. (5.165) simplifies into: 0 b

1 _ rl

_ rv

B v B Gm þ V $ ðbGm 5v m Þ þ V $ Bbhaið1  haiÞ @ vt rm

C C ðv r 5v r ÞC A (5.166)

0 13 X _ X _ 4b ak @sk þ st A5 þ Fkw ¼ brm g  bVP þ V $ k 2

k¼l;v

k¼l;v

where the mixture velocity is given as: vm ¼

_ _ vv Gm ð1  haiÞrl vbl þ hairv c ¼ _ _ rm ð1  haiÞrl þ hairv

(5.167)

and the drift velocity is given as: _

_

vr ¼ vv  vl

(5.168)

Summing up the enthalpy conservation equations (i.e. Eq. 5.140) of the two-fluid model and using the jump interface condition (i.e. Eq. 5.143) lead to: 2 3 2 3



X . . v X4 _ 4b ak r_k hk vc 5 b ak rk hk 5 þ V $ b k vt k¼l;v k¼l;v 2 13 0 2 3 X 6 B _ X _ X

000 00 t 00 C7 t 4b ak s $ vc 5 b a k qk  V $ z 4b ak @qk þ qk A5  V $ k k k¼l;v

k¼l;v

k¼l;v

82 9 0 0 13 2 3 2 1 3 = X

X _ X< _ _ v 4 ak Pk 5 þ b 4 ak @VPk A $ c 4b ak @sk þ st A5: V5 c vk þ b vk 5 þ k : ; vt k¼l;v k¼l;v k¼l;v 2 0 12 3 X X 1 4 Mki @ c þ v k A  Fki $ c v k  Fkw $ c vk 5 þ Ekw 2 k¼l;v k¼l;v (5.169)

Compared to the transfer of energy occurring at the solid walls in nuclear fuel assemblies, it is customary to neglect the following terms in the previous equation (Obry, 1996):

_ P
The volumetric heat sources b ak q000 k . k¼l;v


The work of the forces due to phase change and stresses at the interface P Fki $ c vk . k¼l;v


The work of the forces due to stresses at the solid walls

P k¼l;v

vk . Fkw $ c

Chapter 5
One-/two-phase flow transport and heat transfer 285


The kinetic energy due to the mass variation in case of phase change 0 12 P1 @c vk A . 2 Mki k¼l;v


The viscous and turbulent dissipation terms 82 9 0 2 3 13

<

= P P _ t A5 @ 4 4b ak st $ vc 5 c V $ b . þ : V 5 v a s þ s k k k k k k ; k¼l;v k¼l;v: Assuming uniform pressure distribution between the phases (pressure equilibrium): _

_

_

Pl ¼ Pv ¼ P

(5.170)

and using Eq. (5.144), Eq. (5.169) then simplifies into: 2 3 2 3

X _ . v X4 _ . 5 4b ak rk hk vc 5 b ak rk hk þ V $ b k vt k¼l;v k¼l;v 0 13 9 8 2 3

= X
(5.171)

or 2 b

3 _ rl

_ rv

6 v 6 ðrm hm Þ þ V $ ðbGm hm Þ þ V $ 6bhaið1  haiÞ 4 rm vt

7 7 v r hr 7 5 (5.172)

2

13 0  X
X 6 B _ _ v_ 00 C7 P þ Jm $ VP þ Ekw z  V$ 4b ak @q00k þ qtk A5 þ b vt k¼l;v k¼l;v

where the (static) mixture enthalpy is given as: _ .

hm ¼

_

.

ð1  haiÞrl hl þ hairv hv rm

(5.173)

the enthalpy difference is given as: .

.

hr ¼ hv  hl

(5.174)

and the mixture volumetric flux is given as: vv Jm ¼ ð1  haiÞ vbl þ haic

(5.175)

Although pressure equilibrium between the phases is most often assumed, the two

286

Modelling of Nuclear Reactor Multi-physics

phases are not necessarily in thermal equilibrium (or more exactly in thermodynamic equilibrium). This means that in thermal non-equilibrium, one has: _

_

Tl sTv

(5.176) _

and an extra correlation giving a relationship between Tl

_

and Tv

would also be

required. In the case of thermal equilibrium, one has on the other hand: _

_

Tl ¼ Tv ¼ Tsat

(5.177)

The two phases would thus have the same temperature, equal to the saturation temperature Tsat. In such case, the enthalpy difference as defined by Eq. (5.174) is the latent heat of vaporization. It is interesting at this point to investigate the conditions under which the system of equations in the model presented above can be closed. Listing the number of unknowns and equations (using the same methodology as the one used in Section 5.3.a for the twofluid model), one would notice that there is one unknown too many in order to close the system of equations if both pressure and thermal equilibrium are assumed. In such a case, one additional relationship is needed. Usually, this additional relationship relates the drift velocity to other variables. A common model used to account for different velocities of the two phases is the drift flux model. This model is useful for low pressure and/or low flow rate conditions under steady-state or near-steady-state conditions. The drift flux model can be considered as a topological law. Several drift flux models exist. The most commonly used one is the one suggested by Zuber and Findlay (1965). In this model, the local time-averaged vapour velocity v v is considered as a sum between the two-phase local time-averaged volumetric flux j and a local time-averaged drift velocity of the vapour v vj , i.e. v v ¼ j þ v vj

(5.178)

The two-phase local time-averaged volumetric flux j is given as the sum between the local volumetric fluxes of the liquid and vapour phases as: j ¼ jl þ jv

(5.179)

with the local volumetric flux of the phase k defined as: jk ¼ Xk v k

(5.180)

from which one can demonstrate, using the definition of the time averages given by Eqs. (5.9) and (5.10), that: Z

1 jk ¼ Xk v k ¼ Dt ¼

Dtk 1 Dt Dtk

Z

tþDt 2

tDt 2

1 Xk ðr; tÞv k ðr; tÞdt ¼ Dt

Z

Dtk 2

tþ

Dtk 2

v k ðr; tÞdt

t

(5.181) Dtk 2

tþ

Dtk 2

t

v k ðr; tÞdt ¼ ak v k

Chapter 5
One-/two-phase flow transport and heat transfer 287

where Eq. (5.12) was also used. Applying Eq. (5.181) for the vapour phase gives: jv ¼ av v v

(5.182)

This equation can also be rewritten as: 1

0

jv ¼ av j þ av @v v  jA

(5.183)

Volume-averaging this equation leads to:



hjv i ¼ av j þ av v v  j

(5.184)

Defining the concentration parameter C0 such that:

av j ¼ C0 av hji

(5.185)

and the effective drift velocity Vvj such that: *



av v v  j

+



¼ V vj av

(5.186)

one finds that Eq. (5.184) can be rewritten as:

hjv i ¼ ½C0 hji þ V vj  av

(5.187)

C0 is in the most general case a tensor, whereas Vvj is a vector. Nevertheless, in the case of a mono-directional flow, the tensor C0 has only one non-zero element corresponding to this direction. Furthermore, taking the two-phase average of Eq. (5.181) gives: +

* hjk i ¼

(5.188)

ak v k

If one neglects the spatial and temporal variations of the density on the phase k, one has, according to Eq. (5.124): *

+ ak ðrk v k Þ

vc k ¼



_ ak rk

* 1R 1R ak ðrk v k Þ d 3 r ak v k d 3 r VV VV ¼ z ¼ 1R 1R ak rk d 3 r ak d 3 r VV VV

+ ak v k

ak

(5.189)

Combining Eqs. (5.188) and (5.189), one finds that:

hjk iz ak vc k

(5.190)

288

Modelling of Nuclear Reactor Multi-physics

Taking the two-phase average of Eq. (5.179) thus gives: v l þ hai vc hji ¼ hjl i þ hjv izð1  haiÞ c v ¼ Jm

(5.191)

where the last equality was obtained using Eq. (5.175). Eq. (5.187) can thus be written as: D E hjv i z ½C0 Jm þ V vj  av ¼ ½C0 Jm þ V vj hai

(5.192)

C0 and Vvj are parameters accounting for the spatial distribution of av and j within the flow and of the relative velocity between the phases. Such parameters depend on the flow regime and the pressure. Those parameters are usually given by empirical or semiempirical relationships. As can be seen from Eq. (5.192), the drift flux model gives an additional relationship between Jm, hai and hjv izhai vc v , or equivalently, between Jm, hai and vc v . This additional relationship can be used to close the system of equations of the mixture model with specified drift velocities when pressure equilibrium and thermal equilibrium are assumed. The model above is sometimes referred to as a four-equation model, since it relies on three conservation equations for mixture quantities (mass, momentum and enthalpy, respectively), as well as on a drift flux model.

5.3.c

Homogeneous equilibrium model

The Homogeneous Equilibrium Model (HEM) is a mixture model for which the drift velocity between the phases is identically equal to zero, i.e. _

_

_

_

v r ¼ v v  v l ¼ 05v l ¼ v v ¼ v m

(5.193)

and for which both pressure equilibrium and thermal equilibrium are assumed. Using the equations derived in Section 5.3.b together with Eq. (5.193) leads to:
A mass conservation equation given by: b

v r þ V $ ðbGm Þ ¼ 0 vt m

(5.194)


A momentum conservation equation given by: b

v Gm þ V $ ðbGm 5v m Þ vt _

¼ brm g  bVP þ V $

X

0 13

X _ 4b ak @sk þ st A5 þ Fkw k 2

k¼l;v

(5.195)

k¼l;v


An enthalpy conservation equation given by: v ðr hm Þ þ V $ ðbGm hm Þ vt m 0 2 13  X
X 6 B _ _ v_ 00 C7 z  V$ P þ Jm $ VP þ Ekw 4b ak @q00k þ qtk A5 þ b vt k¼l;v k¼l;v b

(5.196)

Chapter 5
One-/two-phase flow transport and heat transfer 289

where the different quantities are defined as: _

_

rm ¼ ð1  haiÞrl þ hairv

(5.197)

_ _ v v ¼ rm v m Gm ¼ ð1  haiÞrl vbl þ hairv c _ .

hm ¼

_

(5.198)

.

ð1  haiÞrl hl þ hairv hv rm

vv ¼ vm Jm ¼ ð1  haiÞ vbl þ haic

(5.199) (5.200)

where the last equalities in Eqs. (5.198) and (5.200), respectively, were obtained using Eq. (5.193). The system of equations given by Eqs. (5.194)e(5.196), often referred to as a threeequation model, contains as many unknowns as equations and can thus be solved without the need of any extra closure relationship. The HEM is useful for high pressure and high flow rate conditions. The HEM could also be obtained from the mixture models with specified drift velocities derived in Section 5.3.b by setting the parameters in the drift flux model as C0¼I and Vvj ¼ 0.

5.4 Spatial and temporal discretizations of the flow models In this section, the practical implementation of the flow models earlier derived is presented. Due to the variety of flow models available, the spatial and temporal discretization of the space- and time-averaged conservation equations will be presented from a generic viewpoint. In most two-phase flow computer codes used in the nuclear industry, the flow is assumed to be mono-directional along the axial direction z. Considering a fixed control volume V of length Dz having for mid-position the axial elevation z and surface area A(z) (possibly varying along the flow direction), the two-phase volume-average of any quantity c can be rewritten as: hci ¼

¼

1 V

Z

1 Aave Dz

cðr; tÞd 3 r ¼ V

Z

zþDz 2

zDz 2

1 Aave Dz

Z

zþDz 2 zDz 2

Z Aðz0 Þ

cðr; tÞdz0 d 2 r (5.201)

Aðz0 Þfcgdz0

with Aave ¼

1 Dz

Z

zþDz 2

zDz 2

Aðz0 Þdz0

(5.202)

290

Modelling of Nuclear Reactor Multi-physics

and

Z V¼

zþDz 2

zDz 2

Z

dz0 d 2 r ¼

Z

zþDz 2 zDz 2

Aðz0 Þ

Aðz0 Þdz0 ¼ Aave Dz

(5.203)

In Eq. (5.201), {c} represents the two-phase area-average of the quantity c, defined as: fcg ¼

Z

1 AðzÞ

cðr; tÞd 2 r

(5.204)

AðzÞ

Using Eq. (5.201), one then notices that:          v 1 Dz Dz Dz Dz hci ¼ ave fcg z þ A z fcg z  A zþ vz A Dz 2 2 2 2

(5.205)

In such a case and using the two-fluid model as the most general formulation, the space- and time-averaged conservation equations for a mono-directional flow can be transformed into the following ones:
For mass: 0 3 1 1 20 3


_   
_ v4 1 6B Dz Dz 7 _ B cA C cA C b ak rk 5 þ ave 4@bA ak rAk v  @bA ak rAk v 5 ¼ Mki k A zþ k A z vt A Dz 2 2 2

(5.206)


For momentum: 3 2

v4 _ b ak rk c vk 5 vt 9 82 2 0 13 0 13 >  >
_  
_ < 1 Dz Dz = 6 B A c B A c 6 A C7 A C7  4bA ak rAk @ vc þ ave 4bA ak rAk @ vc k 5 v k A5 z þ k 5 v k A5 z  A Dz > 2 2 > ; :

_ zb ak rk g þ

82 9 13 13 0 0 2 >  >  
 _
 _ < 1 Dz Dz = C7 B B 6 6 t;A C7  4b ak @sAk þ st;A 4b ak @sAk þ sk A5 z þ ave k A5 z  A Dz > 2 2 > : ;

2 3     b ak _ _ Dz Dz 7 6 A  PkA z   ave 4Pk z þ 5k þ Fki þ Fkw A Dz 2 2 (5.207)

where k represents the unit vector along the axial direction.

Chapter 5
One-/two-phase flow transport and heat transfer 291


For enthalpy: 3 2

v4 _ . 5 b ak rk hk vt 3 0 1 1
_ . 
_ .    1 6B Dz Dz 7 B C C  @bA ak rAk hAk d v Ak A z þ v Ak A z  þ ave 4@bA ak rAk hAk d 5 A Dz 2 2 20



_ zb ak q000 k 

9 82 13 13 2 0 0 >  >  
 _
 _ = 00 00 1 <6 Dz Dz 00 00 6 B A B A t;A C7 t;A C7  a q þ q a q þ q z þ z  bA bA A 5 A 5 4 @ @ 4 k k k k k k Aave Dz > 2 2 > ; :

0 1 1 3 2 0 13 





  1 6B Dz Dz _ B C C 7 t;A c t;A t A  @bA ak sk $ vc vk  ave 4@bA ak sk $ v Ak A z þ 5 þ 4b ak @sk þ sk A5: V5 c k A z A Dz 2 2 20



3 2 3

    b ak _ _ _ v4 Dz Dz 7 6 c b ak Pk 5 þ ave  PkA z  v k $ 4PkA z þ þ 5k vt 2 2 A Dz 2

0 12 1 v k A  Fki $ c v k  Fkw $ c v k þ Eki þ Ekw þ Mki @ c 2 (5.208)

In the equations above, the superscript A refers to area-averaged quantities defined in a similar manner as the volume-averaged quantities given by Eqs. (5.121)e(5.125). The area-averaged quantities are thus defined as:

_

rAk

8 9 < = R a k rk ak rk d2 r : ; AðzÞ ¼
 ¼ R ak d2 r ak AðzÞ

8 <

_ sAk

(5.209)

9 =

R a s ak sk d 2 r : k k; AðzÞ ¼
 ¼ R ak d 2 r ak AðzÞ

(5.210)

292

Modelling of Nuclear Reactor Multi-physics 8 <

_ PkA

9 =

R a P ak Pk d 2 r : k k; AðzÞ ¼ ¼ R ak d 2 r ak AðzÞ

(5.211)

and 8 < A vc k ¼

:

ak ðrk v k Þ


9 =

R

; ¼

 _ ak rAk

AðzÞ

ak ðrk v k Þ d 2 r R ak rk d 2 r

0

12

B cA C @ vk A . hAk

8 < A ec k ¼

:

_ ¼ ekA



ak ðrk ek Þ

(5.212)

AðzÞ

_

þ

2 9 =


 _ ak rAk

¼

(5.213)

_

rAk

R

;

PkA

ak ðrk ek Þ d 2 r R ak rk d 2 r

AðzÞ

(5.214)

AðzÞ

and 8 > < _ A A b st;A v 0k 5v 0k k z  rk

> : ¼

ak rk v 0k 5v 0k 




9 > = > ; (5.215)

ak 8 > <

_ b

A 00 z rAk ek0 v 0k qt;A k

9 > =  ak rk ek0 v 0k > > : ;
 ¼ ak

(5.216)

It should be emphasized that a similar derivation could be carried out for flows that are not necessarily mono-directional. For the sake of brevity, such a derivation is not presented here, and the interested reader is referred to the existing literature on this subject. The equations presented above (i.e. Eqs. 5.206e5.208) actually correspond to a so-called finite volume discretization scheme, where balance equations were spatially integrated on fixed volumes. As a result of this integration and because of the spatial derivatives, the balance equations contain both volume-averaged and surface-averaged

Chapter 5
One-/two-phase flow transport and heat transfer 293

FIGURE 5.3 Numbering of the control volumes and their boundaries in the finite volume discretization scheme used in modelling flow transport.

quantities. It is interesting to notice that the balance equations written above are, with respect to the spatial discretization, totally equivalent to the nodal balance equations presented in Section 4.2 of Chapter 4 for solving the neutron diffusion equations. Due to the mono-directional assumption of the flow, the discretized system being considered can be sketched in Fig. 5.3. With respect to the control volume i, the volumeaveraged values are denoted with the subscript i, whereas area-averaged values at the boundaries are denoted with the subscripts i  1/2 and i þ 1/2 upstream and downstream, respectively. With the conventions thus defined, the mass, momentum and enthalpy conservation equations can be recast in the following generic form: vui þ vt

h i fiþ12  fi12 Dz

¼ qi

(5.217)

This equation can be considered as a discretized version with respect to space of the following partial differential equation: vu v ðz; tÞ þ f ½uðz; tÞ ¼ qðz; tÞ vt vz

(5.218)

with ui ¼

1 Dz

Z

zþDz 2

zDz 2

uðz0 ; tÞdz0

      Dz Dz and fi12 hf u z  ; t fiþ12 h f u z þ ; t 2 2 qi ¼

1 Dz

Z

zþDz 2

zDz 2

qðz0 ; tÞdz0

(5.219)

(5.220)

(5.221)

since taking the spatial average of Eq. (5.218) leads to Eq. (5.217) with the quantities defined in Eqs. (5.219)e(5.221). In order to find the solution to this system of equations, relationships between ui, fiþ12 and fi12 need to be introduced, i.e. one needs to add extra

294

Modelling of Nuclear Reactor Multi-physics

relationships of the form: fiþ12 ¼ f + ðuiþh ; :::; uihþ1 Þ

(5.222)

with h being a positive integer. This formulation guarantees a conservative discretization scheme, i.e. the fact that the quantity fiþ12 at the interface i þ 1/2 between two neighbouring volumes i and i þ 1 is conserved irrespective of whether the volume i or the volume i þ 1 is considered. The equation above simply states that extra relationships between surface-averaged quantities (denoted with the subscripts i  1/2) and volumeaveraged quantities (denoted with the subscripts i) need to be defined, so that the system of equations to be considered only contains either volume-averaged or surface-averaged quantities. In such a case, the system of equations can be actually solved. The function f+ should be consistent, i.e., should tend to the continuous function f when the mesh size tends to zero. This means that the following property should be fulfilled: f + ðu; :::; uÞ ¼ f ðuÞ

(5.223)

Different possibilities exist. One possibility is to assume that any surface-averaged value is the mean between the two neighbouring node-averaged values. In such a case, one has: fiþ12 ¼ f + ðuiþ1 ; ui Þ ¼

f ðuiþ1 Þ þ f ðui Þ 2

(5.224)

Other possibilities exist, and the interested reader is referred to the literature on that subject for further details (e.g. Eymard et al., 2000). Due to the advection terms in the mass and enthalpy conservation equations, in which the velocities averaged on the surface areas delimiting a node appear, one common approach is to use a spatial discretization scheme according to which the momentum balance equations are solved on volumes shifted by a distance of Dz =2

FIGURE 5.4 Numbering of the control volumes and their boundaries in the finite volume discretization scheme in case of a staggered mesh being used for modelling flow transport.

Chapter 5
One-/two-phase flow transport and heat transfer 295

compared to the discretization on which the mass and enthalpy balance equations are solved, as illustrated in Fig. 5.4. In such a case, the mesh is referred to as a staggered mesh. The advantage of a staggered mesh lies with the fact that the surface-averaged velocities required in the mass and enthalpy conservation equations correspond to the surface-averaged velocities evaluated in the middle of the control volumes in the mesh used for the momentum conservation equation. The surface-averaged velocities in the mesh used for the mass and enthalpy conservation equations are thus replaced by the volume-averaged velocities from the mesh used for the momentum conservation equations. With respect to the velocities and using Eq. (5.222) together with the numbering introduced in Fig. 5.4, this formally reads as: fiþ12 ¼ f + ðuI Þ ¼ f ðuI Þ

(5.225)

Regarding the temporal discretization, different possibilities also exist. Using a firstorder Taylor expansion for the time derivative, the space- and time-discretized version of Eq. (5.217) can be generically written as: h in;nþ1 n fiþ12  fi12 unþ1  u i i þ ¼ qn;nþ1 i Dt Dz

(5.226)

where the subscript n refers to the time-discretized value at time step n of a (spacediscretized) variable. Furthermore, one has for any variable g: g n;nþ1 ¼ qg nþ1 þ ð1  qÞg n with q ˛ ½0; 1

(5.227)

The following time discretization schemes are possible:
Explicit scheme or Euler forward scheme, when q ¼ 0.
Implicit scheme or Euler backward scheme, when q ¼ 1.
Crank-Nicholson scheme, when q ¼ 1=2. In order to investigate the stability of the numerical scheme after the time- and spacediscretizations, it is interesting at this point to turn back to the equivalent partial differential equation given by Eq. (5.218). This equation corresponds to the local conservation equation derived for either mass, momentum or energy (enthalpy), as can be seen from Eqs. (2.73), (2.78) and (2.89) in Chapter 2. If one further assumes, for the purpose of the present discussion, that the flow is incompressible, i.e. V$v¼0

(5.228)

and that one can neglect the source/sink volumetric/surface terms in the local conservation equations, one notices that the local conservation equations for mass and energy (enthalpy) can then be expressed as:

296

Modelling of Nuclear Reactor Multi-physics

vu vu ðz; tÞ þ v ðz; tÞ ¼ 0 vt vz

(5.229)

This equation represents the transport of the quantity u with the flow, or the so-called advection of the quantity u. In case of the Euler forward scheme, and using a first-order Taylor expansion for the spatial derivative, Eq. (5.229) can be written as: unþ1  uni un  uni1 i þv i ¼ 0 if v  0 Dt Dz

(5.230)

un  uni unþ1  uni i þ v iþ1 ¼ 0 if v < 0 Dt Dz

(5.231)

and

whereas in the case of the Euler backward scheme, Eq. (5.229) can be written as: unþ1  uni unþ1  unþ1 i i1 þv i ¼ 0 if v  0 Dt Dz

(5.232)

unþ1  unþ1 unþ1  uni i i þ v iþ1 ¼ 0 if v < 0 Dt Dz

(5.233)

and

In case of positive velocity, the quantity u is transported in the direction of increasing i values, whereas in the case of negative velocity, the quantity u is transported in the direction of decreasing i values. It should be mentioned that it is always necessary to approximate the spatial derivative upwind, i.e. the discretized spatial derivative should contain information from the node upstream from the node being resolved. In such a case, the discretization scheme is compatible with the physical problem being studied. This is referred to as an upwind scheme. Depending on whether the velocity is positive or negative, one thus has two discretization schemes, for both the forward and backward Euler schemes. The discretization scheme defined by Eqs. (5.230) and (5.231) is thus called an upwind explicit scheme or upwind Euler forward scheme, whereas the discretization scheme defined by Eqs. (5.232) and (5.233) is called an upwind implicit scheme or upwind Euler backward scheme. The stability criterion for each of the above two schemes can then be investigated, as explained below. For the purpose of the following discussion, only the case of positive velocities will be considered. The stability of the above discretization schemes can be studied using the Von Neumann method or Fourier method presented in Section 4.4.b of Chapter 4. The solution un+ to the discretized problem can be assumed to contain i round-off errors and thus can be assumed to contain some deviation εni from the exact solution to the discretized problem uni , i.e. un+ ¼ uni þ εni i

(5.234)

Chapter 5
One-/two-phase flow transport and heat transfer 297

Substituting Eq. (5.234) into Eq. (5.230) and Eq. (5.232), for the explicit and implicit and uni fulfil those equations, one schemes, respectively, and using the fact that both un+ i n concludes that the error εi should also fulfil these equations, i.e.
For the explicit scheme: εnþ1  εni εn  εni1 i þv i ¼0 Dt Dz

(5.235)

εnþ1  εni εnþ1  εnþ1 i i1 þv i ¼0 Dt Dz

(5.236)


For the implicit scheme:

Assuming that the computational domain has a length L and N finite volumes of length Dz, the error is developed into a Fourier series as: εni ¼

N X

Anl expðjkl iDzÞ

(5.237)

l¼N

with j¼

pffiffiffiffiffiffiffi 1

(5.238)

and kl ¼

lp NDz

(5.239)

Defining the amplification factor xl such that: Alnþ1 ¼ xl Anl

(5.240)

    maxl xl   1

(5.241)

the stability criterion is given as:

For the upwind explicit scheme, introducing the Fourier expansion given by Eq. (5.237) into Eq. (5.235) leads to: N P l¼N

Anþ1 expðjkl iDzÞ  l

N P l¼N

Anl expðjkl iDzÞ

Dt N P

þv l¼N

Anl expðjkl iDzÞ

(5.242) 

N P l¼N

Dz

Anl exp½jkl ði

 1ÞDz ¼0

or, using Eq. (5.240): N
X l¼N


Anl expðjkl iDzÞ

xl  1 þ

 vDt ½1  expðjkl DzÞ ¼0 Dz

(5.243)

298

Modelling of Nuclear Reactor Multi-physics

from which one concludes that: xl  1 þ

vDt ½1  expðjkl DzÞ ¼ 0; cl Dz

(5.244)

vDt vDt þ expðjkl DzÞ; cl Dz Dz

(5.245)

or xl ¼ 1 

The right-hand side of the above expression actually represents a circle centred at vDt 1  vDt Dz and of radius Dz . By geometric considerations, one could verify that the stability criterion is fulfilled only if: vDt 1 Dz

(5.246)

This thus defines a stability criterion for the upwind explicit scheme, known as the Courant-Friedrichs-Lewy condition or CFL condition, and vDt =Dz is referred to as the Courant number. It thus means that the size of the control volume cannot be too small and/or the time step too large, if one wants the numerical schemes to provide meaningful information, i.e. if one wants the information to be completely ‘advected’ (at the velocity v) within a control volume of size Dz during the time step Dt. Likewise, for the upwind implicit scheme, introducing the Fourier expansion given by Eq. (5.237) into Eq. (5.236) leads to: N P l¼N

Anþ1 expðjkl iDzÞ  l

N P l¼N

Anl expðjkl iDzÞ

Dt N P

þv

l¼N

Anþ1 expðjkl iDzÞ  l

(5.247) N P l¼N

Anþ1 exp½jkl ði  1ÞDz l

Dz

¼0

or, using Eq. (5.240): N
X


Anl expðjkl iDzÞ

l¼N

xl  1 þ

vDt x ½1  expðjkl DzÞ Dz l

 ¼0

(5.248)

from which one concludes that: xl  1 þ

vDt x ½1  expðjkl DzÞ ¼ 0; cl Dz l

(5.249)

or xl ¼

1 ; cl vDt vDt 1þ  expðjkl DzÞ Dz Dz

(5.250)

Fulfilling the stability criterion expressed by Eq. (5.250) is identical in the present case to demonstrating that the norm of the denominator is always larger than unity. The vDt denominator actually represents a circle centred at 1 þ vDt Dz and of radius Dz . By geometric

Chapter 5
One-/two-phase flow transport and heat transfer 299

considerations, one could verify that the stability criterion is always fulfilled, and thus the upwind implicit scheme is stable, whatever size of the control volumes and/or time step one chooses. Finally, for the explicit scheme, one could verify, using a first-order Taylor expansion in time for unþ1 with respect to uni and in space for uni1 with respect to uni , that: i unþ1  uni un  uni1 vu vu i ðzi ; tn Þ þ v ðzi ; tn Þ ¼ 0 þ OðDtÞ þ v i þ OðDzÞ ¼ vt vz Dt Dz

(5.251)

and for the implicit scheme, using a first-order Taylor expansion in time for uni with nþ1 respect to unþ1 and in space for unþ1 , that: i i1 with respect to ui unþ1  uni unþ1  unþ1 vu vu i i1 ðzi ; tnþ1 Þ þ v ðzi ; tnþ1 Þ ¼ 0 þ OðDtÞ þ v i þ OðDzÞ ¼ vt vz Dt Dz

(5.252)

This demonstrates that the error of truncation tends to zero when both the time step and the size of the control volumes tend to zero. This defines a consistent numerical scheme. According to the Lax equivalence theorem, which states that for a well-posed and consistent numerical scheme, stability and convergence are equivalent, one concludes from the above analyses that the implicit scheme is always convergent, whereas the explicit scheme is convergent only when Eq. (5.246) is fulfilled. Although an explicit scheme allows determining a new time estimate of the solution from the previous time estimate in a straightforward manner, it suffers from the necessity to fulfil the CFL condition. The implicit scheme, on the other hand, does not require any condition to be fulfilled. In this respect, implicit schemes have to be preferred to explicit schemes, since they are unconditionally stable. This also means that in case of very fast phenomena, an implicit scheme is necessary in order to prevent severe Courant limitations, i.e. in order to avoid a very small time stepping with an explicit scheme. The implicit scheme nevertheless relies on the source terms in Eq. (5.226), i.e. the terms on the right-hand side of the equation, to be evaluated at the time steps being resolved. This thus requires the inversion of large matrices. The computational effort to find the solution in such a case is much more significant than when using an explicit scheme. When the phenomena being studied are slow, an explicit scheme with a proper size of the control volumes and/or time step might thus be a good enough solution strategy to the problem being investigated. In the previous discussion, quantities (such as mass and energy/enthalpy) advected at the velocity of the fluid were considered. In addition, other phenomena involving for instance pressure waves can propagate through the fluid. Since pressure waves travel at the speed of sound and thus represent very fast phenomena, an implicit scheme is necessary for the discretization of the momentum conservation equation if one wants to avoid severe Courant limitations in those conditions. Due to the variety of the possible combinations of time- and space-discretization, only the basic principles to be followed for studying the stability, convergence and consistency of numerical schemes were introduced.

300

Modelling of Nuclear Reactor Multi-physics

Once the problem has been discretized in space and time, one usually ends up with a large system of algebraic equations, which can be put into a matrix form. The matrix equation is then solved using, e.g., the numerical techniques introduced in Section 4.3 of Chapter 4. Again, due to the variety of possible time- and space-discretization schemes, the interested reader is referred to the existing literature to find out the iteration technique best at solving the system of equations thus obtained. Finally, before closing this section on the spatial and temporal discretization of flow models, the well-posedness character of the two-fluid model is touched upon. A well-posed problem is defined as a problem for which the solution known on a subset of the space-time domain allows determining the solution on the entire space-time domain. Hyperbolic partial differential equations have this character. The local conservation equations are in essence hyperbolic. Nevertheless, the transformations performed on the local conservation equations render the system of equations to be solved non-hyperbolic (Din et al., 2003). The reason for this loss of hyperbolicity lies in the assumption of pressure equilibrium in the derivation of the two-fluid model, as earlier mentioned. When the model being solved becomes ill-posed, numerical instabilities might occur. Several techniques are used at the level of the discretization in order to damp such possible instabilities. The first one is to choose a discretization scheme introducing enough numerical diffusion. The second one is to actually add some artificial terms in the conservation equations of the two-fluid model so that stable solutions are obtained. The issue of well-posedness of the two-fluid model is an intensive subject of research, and the interested reader is thus referred to the existing literature on that subject.

5.5 Modelling of heat conduction in solid structures In this section, heat conduction in solid structures is considered. One of the primary objectives of this section is to derive the model for heat conduction in the fuel elements, although the same model can be used for modelling heat conduction in other solid structures. In most cases, the structures being considered are of cylindrical shapes, in which the radial temperature gradient is far more important than the axial one. For the sake of simplicity of the derivations, emphasis is put on those systems in the following and only the radial heat conduction is considered. The axial dependence is simply treated by dividing the length of the system being modelled into a number of axially averaged regions. In each of these axially averaged regions, the spatial dependence of the variables involved in modelling heat transfer only depends on the distance from the centre of the fuel elements. The heat conduction model is based on the heat conduction equation derived in Chapter 2 (i.e. Eq. (2.88)), and which is recalled hereafter: rðT Þcp ðT Þ

h i vT ðr; tÞ ¼ V $ kðT ÞVT ðr; tÞ þ q000 ðr; tÞ vt

(5.253)

Chapter 5
One-/two-phase flow transport and heat transfer 301

In most applications, the boundary condition associated with this equation is stated using Fourier’s law as: kðT ÞVT ðrB ; tÞ ¼ q00 ðrB ; tÞ

(5.254)

with rB being the position of the boundary. Eq. (5.253) is spatially discretized using a finite volume method, i.e. the equation is integrated on an elementary volume Vi as: Z

rðT Þcp ðT Þ Vi

vT ðr; tÞd 3 r ¼ vt

Z

Z

V $ ½kðT ÞVT ðr; tÞd 3 r þ Vi

q000 ðr; tÞd 3 r

(5.255)

Vi

Using Gauss’ divergence theorem and dividing by the volume Vi, this equation becomes: 1 Vi

Z rðT Þcp ðT Þ Vi

vT 1 ðr; tÞd 3 r ¼ vt Vi

Z kðT ÞVT ðr; tÞ $ Nd 2 r þ Si

1 Vi

Z

q000 ðr; tÞd 3 r

(5.256)

Vi

with Si being the outer surface of the volume Vi and N being the outward normal unit vector to the elementary surface d2r on that surface. If the elementary volume is taken small enough, one can assume that the density, heat capacity and thermal conductivity can be taken at the volume-averaged temperature. If the volume-averaged values are denoted with the subscript i, whereas area-averaged values at the boundaries are denoted with the subscripts i  1/2 and i þ 1/2 in the direction of increasing and decreasing i, respectively, as represented in Fig. 5.5, Eq. (5.256) can be written as: rðTi Þcp ðTi Þ

Si1 vTi Siþ12 z kðTi ÞðVT Þiþ1 þ 2 kðTi ÞðVT Þi1 þ qi 2 2 vt Vi Vi

(5.257)

with Ti ¼

1 Vi

Z T ðr; tÞd 3 r

(5.258)

Vi

FIGURE 5.5 Numbering of the control volumes and their boundaries in the finite volume discretization scheme used in modelling heat conduction.

302

Modelling of Nuclear Reactor Multi-physics

ðVT Þiþ1 ¼ 2

1 Siþ12

Z VT ðr; tÞ $ rd 2 r Siþ1 2

and ðVT Þi1 ¼  2

1 Si12

Z VT ðr; tÞ $ rd 2 r

(5.259)

Si1 2

qi ¼

1 Vi

Z

q000 ðr; tÞd 3 r

(5.260)

Vi

with r being the unit vector oriented in the direction of increasing i indices. In the equations above, the transfer of heat in the axial direction was neglected and thus the surface integrals only refer to integrals taken on surfaces perpendicular to the radial direction, along which a temperature gradient is observed. Denoting by ri the ‘mid-position’ of the region i along the radial direction being modelled, Driþ the distance between the ‘mid-position’ and the boundary to the region i þ 1 and Dri the distance between the ‘mid-position’ and the boundary to the region i  1, as depicted in Fig. 5.6, the gradient terms given by Eq. (5.259) are then discretized as: ðVT Þiþ1;ðiÞ ¼ 2

Tiþ12  Ti Driþ

(5.261)

and

FIGURE 5.6 Geometrical dimensions associated to the region i in the finite volume discretization scheme used in modelling heat conduction.

Chapter 5
One-/two-phase flow transport and heat transfer 303

ðVT Þi1;ðiÞ ¼  2

Ti  Ti12

(5.262)

Dri

The subscript (i) denotes quantities evaluated with respect to the volume i. The approximations of the temperature gradients above rely on the assumption that the temperatures at mid-positions in each of the volumes are equal to the corresponding volume-averaged temperatures. The ‘mid-position’ ri is here defined such that the areas h h i

2 i 2 p ri2  ri  Dri and p ri þ Driþ  ri2 are equal, with ‘mid-position’ referring to a division of the region i into two equi-surface regions. Because of the cylindrical geometry, the distances Dri and Driþ are thus not equal, although they are close to each other. One would also notice that:  

Si12 2 ri þ Driþ 2 ri  Dri ¼ ¼ 2 2 and 2 2 Vi Vi ri þ Driþ  ri  Dri ri þ Driþ  ri  Dri

Siþ12

(5.263)

In order to obtain balance equations only containing volume-averaged quantities, the continuity of the temperatures and heat fluxes at the boundaries i  1/2 and i þ 1/2 will be utilized. Considering first the boundary i  1/2, the temperature gradient can be evaluated in a similar manner as above with respect to the volume i  1 as: ðVT Þi1;ði1Þ ¼ 2

Ti12  Ti1

(5.264)

þ Dri1

The continuity of the heat flux at the boundary i  1/2 imposes that: kðTi1 ÞðVT Þi1;ði1Þ ¼  kðTi ÞðVT Þi1;ðiÞ 2

2

(5.265)

where the minus sign on the right-hand side of the above equation comes from the convention chosen in Eq. (5.259). Using Eqs. (5.262) and (5.264) in Eq. (5.265) leads to: Ti12 ¼

þ kðTi1 ÞDri Ti1 þ kðTi ÞDri1 Ti þ  kðTi1 ÞDri þ kðTi ÞDri1

(5.266)

A similar treatment at the boundary i þ 1/2 would give: Tiþ12 ¼

 kðTi ÞDriþ1 Ti þ kðTiþ1 ÞDriþ Tiþ1  kðTi ÞDriþ1 þ kðTiþ1 ÞDriþ

(5.267)

Using Eqs. (5.261) and (5.262) together with Eqs. (5.266) and (5.267), Eq. (5.257) can be rewritten as: rðTi Þcp ðTi Þ

vTi z ai Ti þ bi Tiþ1 þ ci Ti1 þ qi vt

(5.268)

with ai ¼ kðTi Þ

Siþ12 Vi

Si12 kðTiþ1 Þ kðTi1 Þ þ  kðTi Þ þ þ kðTiþ1 ÞDri Vi kðTi1 ÞDri þ kðTi ÞDri1

 kðTi ÞDriþ1

(5.269)

304

Modelling of Nuclear Reactor Multi-physics

bi ¼ kðTi Þ

ci ¼ kðTi Þ

kðTiþ1 Þ  Vi kðTi ÞDriþ1 þ kðTiþ1 ÞDriþ

(5.270)

kðTi1 Þ þ Vi kðTi1 ÞDri þ kðTi ÞDri1

(5.271)

Siþ12

Si12

It must be noted that at the outer boundaries of the system, the gradient terms in Eq. (5.257) are replaced by:  ðVT Þ1 ¼ 2

1 R 00 q ðrB ; tÞ $ Nd 2 r S12 S1 2

kðT1 Þ

(5.272)

for a boundary located on the left-hand side of the first control volume (i.e. the control volume number 1), and by:  ðVT ÞNþ1 ¼ 2

1 R 00 q ðrB ; tÞ $ Nd 2 r SN þ12 S 1 Nþ

2

kðTN Þ

(5.273)

for a boundary located on the right-hand side of the last control volume (i.e. the control volume number N, in case the system was discretized with N radial nodes). If the lefthand side of the first control volume represents the centre of a fuel rod, one would have ðVT Þ1 ¼ 0. 2 Concerning the gap between the fuel cladding and the fuel pellets in the fuel elements, the surface heat flux at an intermediate gap position is often given by q00 ¼ hg ðTfo  Tci Þ

(5.274)

where hg is the gap conductance, Tfo is the outer temperature of the fuel pellet and Tci is the inner temperature of the cladding (Todreas and Kazimi, 1993). The heat transfer mechanisms occurring in the gap involve transfer by conduction and by radiation. Due to the smallness of the gap, the heat transfer can be alternatively modelled by setting a thermal conductivity kg for the gap equal to: kg ¼ hg d

(5.275)

where d is the thickness of the gap and with a density rg equal to zero. It thus means that the heat transfer through the gap can be equivalently modelled as a heat transfer by conduction. Consequently, the heat transfer through the gap can also be modelled using the same balance equation as the one used for modelling the heat transfer in the fuel pellet and in the fuel cladding, i.e. using Eq. (5.268) without heat source. Regarding the temporal discretization, different possibilities exist. Using a first-order Taylor expansion for the time derivative, the time- and space-discretized version of Eq. (5.253) can be generically written as:   T nþ1  Tin

n;nþ1 n;nþ1 z an;nþ1 r Tin;nþ1 cp Tin;nþ1 i Tin;nþ1 þ bn;nþ1 Tiþ1 þ cin;nþ1 Ti1 þ qn;nþ1 i i i Dt

(5.276)

Chapter 5
One-/two-phase flow transport and heat transfer 305

where the subscript n refers to the time-discretized value at time step n of a (spacediscretized) variable. Furthermore, one has for any variable g: g n;nþ1 ¼ qg nþ1 þ ð1  qÞg n with q ˛ ½0; 1

(5.277)

As for fluid transport modelling, the following time discretization schemes are possible:
Explicit scheme or Euler forward scheme, when q ¼ 0.
Implicit scheme or Euler backward scheme, when q ¼ 1.
Crank-Nicholson scheme, when q ¼ 1=2. In order to investigate the stability of the numerical scheme after the time- and spacediscretization, the volumetric heat source is neglected, and one assumes temperatureindependent thermal properties. Eq. (5.276) can then be recast into: Tinþ1  Tin ðTiþ1  Ti Þn;nþ1 ðTi  Ti1 Þn;nþ1  a za Dt ðDrÞ2 ðDrÞ2

(5.278)

with a¼

k 4rcp

(5.279)

Eqs. (5.278) and (5.279) were derived assuming that ri [Driþ and ri [Dri and that the node spacing was uniform, i.e. Driþ ¼ Dr ¼ Dri ; ci. In such conditions, the problem being modelled is actually equivalent to a homogeneous one-dimensional slab and Eq. (5.278) is the time- and space-discretized version of the following differential equation: vT ðr; tÞ ¼ aV2 T ðr; tÞ vt

(5.280)

In case of the Euler forward scheme, Eq. (5.278) can be written as: n Tn  Tn Tinþ1  Tin T n  Ti1 z a iþ1 2 i  a i 2 Dt ðDrÞ ðDrÞ

(5.281)

In case of the Euler backward scheme, Eq. (5.278) can be written as: nþ1 T nþ1  Tinþ1 Tinþ1  Tin T nþ1  Ti1 z a iþ1 a i 2 2 Dt ðDrÞ ðDrÞ

(5.282)

The stability of the above discretization schemes can be studied using the Von Neumann method or Fourier method presented in Section 4.4.b of Chapter 4. The solution un+ to the discretized problem can be assumed to contain round-off errors, and i thus can be assumed to contain some deviation εni from the exact solution to the discretized problem uni , i.e. un+ ¼ uni þ εni i

(5.283)

Substituting Eq. (5.283) into Eq. (5.281) and Eq. (5.282), for the explicit and implicit

306

Modelling of Nuclear Reactor Multi-physics

schemes, respectively, and using the fact that both un+ and uni fulfil those equations, one i n concludes that the error εi should also fulfil these equations, i.e.
For the explicit scheme: εn  εn εnþ1  εni εn  εni1 i z a iþ1 2 i  a i Dt ðDrÞ ðDrÞ2

(5.284)

εnþ1  εnþ1 εnþ1  εni εnþ1  εnþ1 i i i1 z a iþ1 a i 2 Dt ðDrÞ ðDrÞ2

(5.285)


For the implicit scheme:

Assuming that the computational domain has a length L discretized with N finite volumes of length Dz, the error is developed into a Fourier series as: εni ¼

N X

Anl expðjkl iDzÞ

(5.286)

l¼N

with j¼

pffiffiffiffiffiffiffi 1

(5.287)

and kl ¼

lp NDz

(5.288)

Defining the amplification factor xl such that: Anþ1 ¼ xl Anl l

(5.289)

  maxl xl   1

(5.290)

the stability criterion is given as:

For the explicit scheme, introducing the Fourier expansion given by Eq. (5.286) into Eq. (5.284) leads to: xl ¼ 1 þ 2a

Dt ðDrÞ2

½cosðkl DzÞ  1; cl

(5.291)

The maximum value of xl is 1, whereas its minimum value is given by . 1  4aDt ðDrÞ2 . Fulfilling the stability criterion expressed by Eq. (5.290) thus only imposes a condition on the lower bound, from which one concludes that:

Chapter 5
One-/two-phase flow transport and heat transfer 307

a

Dt ðDrÞ2



1 2

(5.292)

For the implicit scheme, introducing the Fourier expansion given by Eq. (5.286) into Eq. (5.285) leads to: xl ¼

1

1  2a

Dt ðDrÞ

; cl ½cosðk DzÞ  1 l 2

(5.293)

The maximum value of xl is 1, whereas its minimum value is given by i1 h . 1 þ 4aDt ðDrÞ2 , which is positive. Therefore, the stability criterion expressed by Eq. (5.290) is always fulfilled. Finally, for the explicit scheme, one could verify, using a first-order Taylor expansion n in time for Tinþ1 with respect to Tin and a second-order Taylor expansion in space for Tiþ1 n n and Ti1 with respect to Ti , that: " # n n Tiþ1  Tin Tinþ1  Tin Tin  Ti1 þ OðDtÞ  a a þ OðDrÞ Dt ðDrÞ2 ðDrÞ2

(5.294)

vT ðri ; tn Þ  aV2 T ðri ; tn Þ ¼ 0 ¼ vt

and for the implicit scheme, using a first-order Taylor expansion in time for Tin with nþ1 nþ1 respect to Tinþ1 and a second-order Taylor expansion in space for Tiþ1 and Ti1 with respect to Tinþ1 , that:

" # nþ1 nþ1 Tiþ1  Tinþ1 Tinþ1  Tin Tinþ1  Ti1 þ OðDtÞ  a a þ OðDrÞ Dt ðDrÞ2 ðDrÞ2

(5.295)

vT ðri ; tnþ1 Þ  aV2 T ðri ; tnþ1 Þ ¼ 0 ¼ vt

This demonstrates that the error of truncation tends to zero when both the time step and the size of the control volumes tend to zero. This defines a consistent numerical scheme. According to the Lax equivalence theorem, which states that for a well-posed and consistent numerical scheme, stability and convergence are equivalent, one concludes from the above analyses that the implicit scheme is always convergent, whereas the explicit scheme is convergent only when Eq. (5.292) is fulfilled.

5.6 Conclusions In this chapter, the macroscopic modelling of nuclear thermal-hydraulic phenomena was presented. Concerning the modelling of the flow, a set of time-averaged and volume-averaged balance equations was derived for expressing the conservation of mass, momentum

308

Modelling of Nuclear Reactor Multi-physics

and energy/enthalpy. The length of the time-averaging windows and the size of the control volume are chosen so that the number of unknowns is compatible with the physics of the phenomena being studied while guaranteeing a stable, consistent and convergent numerical solution. Due to the size of nuclear reactors, the nuclear fuel assemblies are usually represented by a set of axially stacked control volumes, having for radial size the size of a fuel assembly and for axial size ca. 1/20 of the axial length of a fuel assembly. It should be noted that some nuclear thermal-hydraulic codes, called subchannel codes, have the ability to use as control volumes, an array of fuel pins with their surrounding coolant or even single fuel pins with their respective coolants. The macroscopic approach used results in a filtering of high-frequency and smallscale phenomena and in the introduction in the balance equations of interfacial terms and wall terms. Such terms can only be known experimentally. It also means that most of the macroscopic flow models heavily rely on experimentally derived correlations necessary to close the system of equations. Although the balance equations were derived in the arbitrary case of the possible simultaneous presence of both a liquid and vapour phase, such balance equations can be readily used when only one phase is present as well. It should also be mentioned that some nuclear thermal-hydraulic multi-field codes can handle more than two ‘fields’. In the case of, e.g., an annular flow regime, three fields could be modelled: the liquid film, the vapour core and the liquid droplets dispersed in the vapour core. For each field, a set of balance equations expressing the conservation of mass, momentum and energy/ enthalpy needs to be formulated. The modelling of heat conduction in the fuel elements is, in appearance, considerably simpler than the flow modelling. The difficulty in the modelling of heat conduction lies with obtaining the parameters appearing in the heat conduction equation, such as the densities, heat capacities and thermal conductivities. Such parameters heavily rely on the fuel burn-up and utilization (Graves, 1979). The fuel pellets experience thermal distortion, pellet cracking, fuel densification and swelling. Fission gases are also released during irradiation. In addition, contact points between the fuel pellet and the cladding due to gap closure during operation and irradiation create fuel pellet/cladding interaction stresses. Finally, the cladding itself is exposed to thermal and pressure stresses. All these phenomena, having an impact on the parameters appearing in the heat conduction equation, make the determination of those parameters far from trivial. It should be noted that such phenomena are accounted for in fuel performance codes, which are codes specifically developed for modelling the behaviour of the fuel elements during their utilization.

References Analytis, G.T., 2003. Lectures on Nuclear Engineering. Lecture Notes. Chalmers University of Technology, Gothenburg, Sweden.

Chapter 5
One-/two-phase flow transport and heat transfer 309

Berry, R.A., 2003. Ensemble Averaged Conservation Equations for Multiphase, Multi-Component, and Multi-Material Flows. Technical Report INEEL/EXT-03-01011. Idaho National Engineering and Environmental Laboratory, USA. Delhaye, J.M., 2008. Thermohydraulique des re´acteurs. EDP Sciences, Les Ulis, France. Dinh, T.N., Nourgaliev, R.R., Theofanous, T.G., 2003. Understanding the ill-posed two-fluid model. In: Proceedings of the 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics. NURETH-10, 5 e 9 October 2003, Seoul, Korea. Korean Nuclear Society, Taejon, Korea. Drew, D.A., Passman, S.L., 1999. Theory of Multicomponent Fluids. Springer-Verlag, New York, USA. Eymard, R., Galloue¨t, T., Herbin, R., 2000. Finite volume methods. In: Ciarlet, P.G., Lions, J.L. (Eds.), Techniques of Scientific Computing, Part III, Handbook of Numerical Analysis, VII. North-Holland, Amsterdam, the Netherlands. Graves, H.W., 1979. Nuclear Fuel Management. John Wiley & Sons, Inc., New York, USA. Hewitt, G.F., Roberts, D.N., 1969. Studies of Two-phase Flow Patterns by Simultaneous X-Ray and Flash Photography. Technical report AERE-M2159. United Kingdom Atomic Energy Authority, UK. Ishii, M., Hibiki, T., 2011. Thermo-fluid Dynamics of Two-phase Flow, second ed. Springer Science þ Business Media, LCC, New York, USA. Obry, P., 1996. Ele´ments sur les e´coulements diphasiques. Lecture notes, Institut National des Sciences et Techniques Nucle´aires (INSTN), Saclay, France. Taitel, Y., Dukler, A.E., 1977. Flow regime transitions for vertical upward gas-liquid low: a preliminary approach. In: Proceedings of the 70th Annual Meeting of the AIChE, November 13 e 17, 1977, New York, USA. American Institute of Chemical Engineers, USA. Taitel, Y., Bornea, D., Dukler, A.E., 1980. Modelling flow pattern transitions for steady upward gas-liquid flow in vertical tubes. AIChE Journal 26 (3), 345e354. Todreas, N.E., Kazimi, M.S., 1993. Nuclear Systems I: Thermal Hydraulic Fundamentals. Taylor & Francis, Levittown, USA. Zuber, N., Findlay, J.A., 1965. Average volumetric concentration in two-phase flow systems. Journal of Heat Transfer 87 (4), 453e468.