Difference between revisions of "Kharon theory manual"

Revision as of 11:24, 1 April 2019

1 Conservation of Mass
2 Conservation of Momentum
3 Conservation of Total Enthalpy
4 Face Velocity Equation
5 Phase-Coupled Face Velocity Equation
6 Pressure & Volume Fraction Correction
7 Solution Algorithm
8 Momentum Transfer
- 8.1 Pressure loss
- 8.2 Interface momentum transfer
9 Heat and Mass Transfer
10 References

Conservation of Mass

Figure 1: Conservation of mass for an interior cell

Steady one-dimensional conservation of mass (of phase q) is given through the following relation

$\nabla\cdot\left(\varepsilon \alpha_q \rho_q u_q\right) = \gamma_{pq} + s_{\text{mass},q}$

(1)

Equation 1 is discretized for an interior cell, shown in Figure 1, by volume integration over cell P.

$\Leftrightarrow \iiint\limits_V \left[\nabla\cdot\left(\varepsilon \alpha_q \rho_q u_q\right)\right]\, \text{d}V = \iiint\limits_V \left(\gamma_{pq} + s_{\text{mass},q}\right)\, \text{d}V$

(2)

The divergence term on the left-hand side of the equation is transformed into a surface integral over the cell faces using the divergence theorem.

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q u_q\right)\cdot n\right]\, \text{d}A = \iiint\limits_V \left(\gamma_{pq} + s_{\text{mass},q}\right)\, \text{d}V$

(3)

Throughout this document upwind discretization is used to evaluate the face mass flow rates, even though any discretization scheme could be chosen. The final discretized form reduces to

$\Leftrightarrow \left(\varepsilon \alpha_q \rho_q u_q A\right)_\text{e} - \left(\varepsilon \alpha_q \rho_q u_q A\right)_\text{w} = F_{\text{e},q} - F_{\text{w},q} = \Gamma_{pq} + S_{\text{mass},q}$

(4)

where:

$\alpha_q$	volume fraction of phase q [-],
$\gamma_{pq}$	volumetric mass transfer rate from phase p to q [kg/m³s],
$\Gamma_{pq}$	mass transfer rate from phase p to q [kg/s],
$\varepsilon$	porosity, fluid fraction of cell volume [-],
$\rho_q$	density of phase q [kg/m³],
$A$	face area [m²],
$F$	face mass flow rate [kg/s],
$n$	face normal pointing out of the cell (1 for east face, -1 for west face),
$s_{\text{mass},q}$	volumetric mass source to phase q [kg/m³s],
$S_{\text{mass},q}$	mass source to phase q [kg/s],
$u_q$	velocity of phase q [m/s],
subscript e	east face (the face in the negative x-direction),
subscript w	west face (the face in the positive x-direction),
subscript q	current phase for which the equation is written (1 = primary, 2 = secondary),
subscript p	the other phase (p = 2 for q = 1 and p = 1 for q = 2),
subscript pq	indicates exchange between phases (e.g. $\Gamma_{pq}$ mass transfer from phase p to q).

Conservation of Momentum

Figure 2: Conservation of momentum for an interior cell

Steady one-dimensional conservation of momentum (of phase q) is solved from the following relation

$\nabla\cdot\left(\varepsilon \alpha_q \rho_q u_q u_q\right) = \nabla\cdot\left(\varepsilon \alpha_q \mu_{\text{eff},q} \nabla u_q\right)$

(5)

+ \left[\text{max}\left(\gamma_{pq},0\right)u_p - \text{max}\left(-\gamma_{pq},0\right)u_q\right] - \varepsilon \alpha_q \nabla p + \varepsilon \alpha_q \rho_q g

- \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C \left| u_q \right| u_q - k\left(u_q - u_p\right) + s_{\text{mom},q}

Equation 5 is discretized for an interior cell, shown in Figure 2, by volume integration over cell P. The diffusion term has been moved to the left-hand side, since it will be developed together with the divergence of momentum.

$\Leftrightarrow \iiint\limits_V \left[\nabla\cdot\left(\varepsilon \alpha_q \rho_q u_q u_q\right) - \nabla\cdot\left(\varepsilon \alpha_q \mu_{\text{eff},q} \nabla u_q\right)\right]\, \text{d}V = \iiint\limits_V \left[\text{max}\left(\gamma_{pq},0\right)u_p - \text{max}\left(-\gamma_{pq},0\right)u_q\right]\, \text{d}V$

(6)

+ \iiint\limits_V \left(-\varepsilon \alpha_q \nabla p + \varepsilon \alpha_q \rho_q g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C \left| u_q \right| u_q - k\left(u_q - u_p\right) + s_{\text{mom},q}\right)\, \text{d}V

As with the conservation of mass equation above, the divergence terms on the left-hand side of the equation are transformed into surface integrals over the cell faces using the divergence theorem.

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q u_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x}\right)\cdot n\right]\, \text{d}A = \left[\text{max}\left(\Gamma_{pq},0\right)u_p - \text{max}\left(-\Gamma_{pq},0\right)u_q\right]$

(7)

-\varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q - K\left(u_q - u_p\right) + S_{\text{mom},q}

Developing the left-hand side terms first

$\iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q u_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(8)

= \left(\varepsilon \alpha_q \rho_q u_q u_q A\right)_{\text{e}} - \left(\varepsilon \alpha_q \rho_q u_q u_q A\right)_{\text{w}} - \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x} A\right)_{\text{e}} + \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x} A\right)_{\text{w}}

= F_{\text{e},q} u_{\text{e},q} + F_{\text{w},q} u_{\text{w},q} - \left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{e}} \left(u_{\text{e},q} - u_{\text{P},q}\right) + \left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{w}} \left(u_{\text{P},q} - u_{\text{w},q}\right)

Adopting upwind discretization for the momentum fluxes on the cell faces

$F_{\text{e},q} u_{\text{e},q} = \text{max}\left(F_{\text{e},q},0\right)u_{\text{P},q} - \text{max}\left(-F_{\text{e},q},0\right)u_{\text{E},q}$

(9)

$F_{\text{w},q} u_{\text{w},q} = \text{max}\left(F_{\text{w},q},0\right)u_{\text{W},q} - \text{max}\left(-F_{\text{w},q},0\right)u_{\text{P},q}$

(10)

and rearranging terms

$\Leftrightarrow\iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q u_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(11)

=\left[\text{max}\left(F_{\text{e},q},0\right) + \text{max}\left(-F_{\text{w},q},0\right) + \left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{e}} + \left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{w}}\right]u_{\text{P},q}

-\left[\text{max}\left(-F_{\text{e},q},0\right) + \left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{e}}\right]u_{\text{E},q}

-\left[\text{max}\left(F_{\text{w},q},0\right) + \left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{w}}\right]u_{\text{W},q}

Denoting $\left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{e}} = D_{\text{e},q}$ and $\left(\varepsilon \alpha_q A \frac{\mu_{\text{eff},q}}{\delta x}\right)_{\text{w}} = D_{\text{w},q}$ and modifying the central coefficient slightly, since $\text{max}\left(F_{\text{e},q},0\right) = F_{\text{e},q} + \text{max}\left(-F_{\text{e},q},0\right)$ and $\text{max}\left(-F_{\text{w},q},0\right) = -F_{\text{w},q} + \text{max}\left(F_{\text{w},q},0\right)$ , Equation 11 reduces to

$\Leftrightarrow\iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q u_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(12)

= \left[\text{max}\left(-F_{\text{e},q},0\right) + D_{\text{e},q} + \text{max}\left(F_{\text{w},q},0\right) + D_{\text{w},q} + \left(F_{\text{e},q} - F_{\text{w},q}\right)\right]u_{\text{P},q}

-\left[\text{max}\left(-F_{\text{e},q},0\right) + D_{\text{e},q}\right]u_{\text{E},q} - \left[\text{max}\left(F_{\text{w},q},0\right) + D_{\text{w},q}\right]u_{\text{W},q}

Identifying the coefficients of neighboring velocities

$a_{\text{e},q} = \text{max}\left(-F_{\text{e},q},0\right) + D_{\text{e},q}$

(13)

$a_{\text{w},q} = \text{max}\left(F_{\text{w},q},0\right) + D_{\text{w},q}$

(14)

leads to a general form

$\Leftrightarrow\iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q u_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \mu_{\text{eff},q} \frac{\Delta u_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(15)

= \left[a_{\text{e},q} + a_{\text{w},q} + \left(F_{\text{e},q} - F_{\text{w},q}\right)\right]u_{\text{P},q} - a_{\text{e},q}u_{\text{E},q} - a_{\text{w},q}u_{\text{W},q}

Substituting the results of Equation 15 into Equation 7 and continuing with the development

$\Rightarrow\left[a_{\text{e},q} + a_{\text{w},q} + \left(F_{\text{e},q} - F_{\text{w},q}\right)\right]u_{\text{P},q} - a_{\text{e},q}u_{\text{E},q} - a_{\text{w},q}u_{\text{W},q}$

(16)

= \left[\text{max}\left(\Gamma_{pq},0\right)u_p - \text{max}\left(-\Gamma_{pq},0\right)u_q\right] - \varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g

- \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q - K\left(u_q - u_p\right) + S_{\text{mom},q}

At this point, the conservation of mass equation times cell velocity, $u_{\text{P},q}$ , is subtracted from Equation 16

$\Rightarrow\left[a_{\text{e},q} + a_{\text{w},q} + \left(F_{\text{e},q} - F_{\text{w},q}\right) - \left(F_{\text{e},q} - F_{\text{w},q}\right)\right]u_{\text{P},q} - a_{\text{e},q}u_{\text{E},q} - a_{\text{w},q}u_{\text{W},q}$

(17)

= \left[\text{max}\left(\Gamma_{pq},0\right)u_p - \text{max}\left(-\Gamma_{pq},0\right)u_q\right] - \left(\Gamma_{pq} + S_{\text{mass},q}\right)u_{\text{P},q}

-\varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q - K\left(u_q - u_p\right) + S_{\text{mom},q}

Splitting the mass source term

$\Leftrightarrow\left[a_{\text{e},q} + a_{\text{w},q}\right]u_{\text{P},q} - a_{\text{e},q}u_{\text{E},q} - a_{\text{w},q}u_{\text{W},q}$

(18)

= \left[\text{max}\left(\Gamma_{pq},0\right)u_{\text{P},p} - \text{max}\left(-\Gamma_{pq},0\right)u_{\text{P},q}\right]

-\left[\text{max}\left(\Gamma_{pq},0\right) - \text{max}\left(-\Gamma_{pq},0\right) + \text{max}\left(S_{\text{mass},q},0\right) - \text{max}\left(-S_{\text{mass},q},0\right)\right]u_{\text{P},q}

-\varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q - K\left(u_q - u_p\right) + S_{\text{mom},q}

$\Leftrightarrow\left[a_{\text{e},q} + a_{\text{w},q} + \text{max}\left(\Gamma_{pq},0\right) + \text{max}\left(S_{\text{mass},q},0\right) + \text{max}\left(-\Gamma_{pq},0\right)\right]u_{\text{P},q} - a_{\text{e},q}u_{\text{E},q} - a_{\text{w},q}u_{\text{W},q}$

(19)

= \left[\text{max}\left(-\Gamma_{pq},0\right) + \text{max}\left(-S_{\text{mass},q},0\right)\right]u_q + \text{max}\left(\Gamma_{pq},0\right)u_p

-\varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q - K\left(u_q - u_p\right) + S_{\text{mom},q}

Since the right-hand side contains only references to the central cell, subscripts P have been dropped from right-hand side to simplify the notation. Note that the diagonal element of the coefficient matrix contains the split terms, $\text{max}\left(\Gamma_{pq},0\right) + \text{max}\left(S_{\text{mass},q},0\right)$ , resulting from the subtraction of the mass equation times cell velocity, while the counterparts, $\left[\text{max}\left(-\Gamma_{pq},0\right) + \text{max}\left(-S_{\text{mass},q},0\right)\right]u_q$ , remain on the right-hand side. These are not considered as momentum source terms, as such, since they are only artifacts of the derivation procedure.

In preparation for casting Equation 19 into a vector form over the phases, a phase coupling matrix, $\mathbf{K}$ , is formed

$\begin{cases} \mathbf{K} = \begin{bmatrix} K & -K \\ -K -\Gamma_{12} & K + \Gamma_{12} \end{bmatrix} \qquad \Gamma_{12} \ge 0 \\ \\ \mathbf{K} = \begin{bmatrix} K - \Gamma_{12} & -K + \Gamma_{12} \\ -K & K \end{bmatrix} \qquad \Gamma_{12} < 0 \end{cases}$

(20)

Finally, denoting the central coefficient by $a_P = a_{\text{e},q} + a_{\text{w},q} + \text{max}\left(S_{\text{mass},q},0\right)$ and casting Equation 19 into a vector form (over the phases) the final form, in which the equations are actually solved, is obtained. For simplicity’s sake, the equation is written only for the first 3 cells (c1, c2 and c3).

$\begin{bmatrix} a_{\text{P,c1},1} + \mathbf{K}_{\text{11,c1}} & \mathbf{K}_{\text{12,c1}} & -a_{\text{e,c1},1} & 0 & \cdots & 0 \\ \mathbf{K}_{\text{21,c1}} & a_{\text{P,c1},2} + \mathbf{K}_{\text{22,c1}} & 0 & -a_{\text{e,c1},2} & \ddots & \vdots \\ -a_{\text{w,c2},1} & 0 & a_{\text{P,c2},1} + \mathbf{K}_{\text{11,c2}} & \mathbf{K}_{\text{12,c2}} & -a_{\text{e,c2},1} & 0 \\ 0 & -a_{\text{w,c2},2} & \mathbf{K}_{\text{21,c2}} & a_{\text{P,c2},2} + \mathbf{K}_{\text{22,c2}} & 0 & -a_{\text{e,c2},2} \\ \vdots & \ddots & -a_{\text{w,c3},1} & 0 & a_{\text{P,c3},1} + \mathbf{K}_{\text{11,c3}} & \mathbf{K}_{\text{12,c3}} \\ 0 & \cdots & 0 & -a_{\text{w,c3},2} & \mathbf{K}_{\text{21,c3}} & a_{\text{P,c3},2} + \mathbf{K}_{\text{22,c3}} \end{bmatrix} \begin{pmatrix} u_{\text{c1,1}}\\ u_{\text{c1,2}}\\ u_{\text{c2,1}}\\ u_{\text{c2,2}}\\ u_{\text{c3,1}}\\ u_{\text{c3,2}} \end{pmatrix} = \begin{pmatrix} S_{\text{mom,c1,1}}\\ S_{\text{mom,c1,2}}\\ S_{\text{mom,c2,1}}\\ S_{\text{mom,c2,2}}\\ S_{\text{mom,c3,1}}\\ S_{\text{mom,c3,2}} \end{pmatrix}$

(21)

where:

S_{\text{mom},\text{c1},q} = \Big[ \text{max}\left(-S_{\text{mass},q},0\right)u_q -\varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q + S_{\text{mom},q} \Big]_{\text{c1}}

and

$\alpha_q$		volume fraction of phase q [-],
$\gamma_{pq}$		volumetric mass transfer rate from phase p to q [kg/m³s],
$\Gamma_{pq}$		mass transfer rate from phase p to q [kg/s],
$\delta x$		distance between two adjacent cell centers [m],
$\varepsilon$		porosity, fluid fraction of cell volume [-],
$\mu_{\text{eff},q}$		effective viscosity of phase q [Ns/m²],
$\rho_q$		density of phase q [kg/m³],
$A$		face area [m²],
$C$		inertial loss coefficient [1/m],
$F$		face mass flow rate [kg/s],
$g$		acceleration of gravity [m/s²],
$k$		volumetric interphase drag coefficient [kg/m³s],
$K$		interphase drag coefficient [kg/s],
$n$		face normal pointing out of the cell (1 for east face, -1 for west face),
$p$		pressure [Pa],
$s_{\text{mass},q}$		volumetric mass source to phase q [kg/m³s],
$S_{\text{mass},q}$		mass source to phase q [kg/s],
$s_{\text{mom},q}$		volumetric momentum source to phase q [N/m³],
$S_{\text{mom},q}$		momentum source to phase q [N],
$u_q$		velocity of phase q [m/s],
$V$		cell volume [m³],
subscript c1		refers to first cell,
subscript c2		refers to second cell,
subscript c3		refers to third cell,
subscript e		east face (the face in the negative x-direction),
subscript E		east neighbor (the cell in the negative x-direction),
subscript w		west face (the face in the positive x-direction),
subscript W		west neighbor (the cell in the positive x-direction),
subscript P		central cell (the cell for which the equation is written),
subscript q		current phase for which the equation is written (1 = primary, 2 = secondary),
subscript p		the other phase (p = 2 for q = 1 and p = 1 for q = 2),
subscript pq		indicates exchange between phases (e.g. $\Gamma_{pq}$ mass transfer from phase p to q).

Conservation of Total Enthalpy

Figure 3: Conservation of total enthalpy for an interior cell

Conservation of total enthalpy (of phase q) can be written as

$\frac{\partial \left(\varepsilon \alpha_q \rho_q h_{0,q}\right)}{\partial t} + \nabla\cdot\left(\varepsilon \alpha_q \rho_q h_{0,q} u_q\right)$

(22)

= \nabla\cdot\left(\varepsilon \alpha_q \lambda_{\text{eff},q} \nabla T_q \right) + \left[\text{max}\left(\gamma_{pq},0\right) h_{0,\text{cre}} - \text{max}\left(-\gamma_{pq},0\right) h_{0,\text{des}}\right]

- \nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \varepsilon \alpha_q \frac{\partial p}{\partial t} + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q}

Disregarding the transient terms and moving the diffusion term to the left-hand side leads to

$\Rightarrow\nabla\cdot\left(\varepsilon \alpha_q \rho_q h_{0,q} u_q\right) - \nabla\cdot\left(\varepsilon \alpha_q \lambda_{\text{eff},q} \nabla T_q \right)$

(23)

= \left[\text{max}\left(\gamma_{pq},0\right) h_{0,\text{cre}} - \text{max}\left(-\gamma_{pq},0\right) h_{0,\text{des}}\right] - \nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q}

Total enthalpy is defined as the sum of specific enthalpy, kinetic energy and potential energy

$h_0 = h + \tfrac{1}{2} \left| u_q \right|^2 - \mathbf{g} \cdot \mathbf{r}$

(24)

where $\mathbf{g} \cdot \mathbf{r}$ is the dot product of gravity vector, $\mathbf{g}$ , and location vector, $\mathbf{r}$ , a vector pointing from the reference level to the evaluated location. The total enthalpies related to phase change can be defined as

$\begin{align} & h_{0,\text{cre}} = h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\\ & h_{0,\text{des}} = h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}} \end{align}$

(25)

In other words, when mass is transferred to phase q (“created”), it brings with it the saturation enthalpy of phase q, $h_{\text{P,sat,}q}$ , and the kinetic energy of phase p, $\tfrac{1}{2} \left| u_p \right|^2_{\text{P}}$ . When mass is transferred from phase q (“destroyed”), it leaves with the enthalpy and kinetic energy of phase q, $h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}}$ . It is chosen that the potential energy associated with phase change is evaluated at the center of the cell, $\left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}$ .

Substituting Equation 24 and Equation 25 to Equation 23 and rearranging terms

$\Rightarrow\nabla\cdot\left(\varepsilon \alpha_q \rho_q h_q u_q\right) - \nabla\cdot\left(\varepsilon \alpha_q \lambda_{\text{eff},q} \nabla T_q \right)$

(26)

= \left[\text{max}\left(\gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}} \right)\right]

-\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right) +\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right) -\nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q}

At this point, an approximation is introduced: the diffusion term is expressed relative to the gradient of enthalpy, instead of temperature, as follows. This approximation is considered to be sufficiently accurate for small enthalpy differences.

$\Rightarrow\nabla\cdot\left(\varepsilon \alpha_q \rho_q h_q u_q\right) - \nabla\cdot\left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \nabla h_q \right)$

(27)

= \left[\text{max}\left(\gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right]

-\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right) +\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right) -\nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q}

Equation 27 is discretized for an interior cell, shown in Figure 3, by volume integration over cell P.

$\Leftrightarrow \iiint\limits_V \left[\nabla\cdot\left(\varepsilon \alpha_q \rho_q h_q u_q\right) - \nabla\cdot\left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \nabla h_q \right)\right]\, \text{d}V$

(28)

= \iiint\limits_V \Big\{ \left[\text{max}\left(\gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right]

-\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right) +\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right) -\nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q} \Big\}\, \text{d}V

As with the conservation of mass and momentum above, the divergence terms on the left-hand side of the equation are transformed into surface integrals over the cell faces using the divergence theorem.

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q h_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(29)

= \iiint\limits_V \Big\{ \left[\text{max}\left(\gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right]

-\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right) +\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right) -\nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q} \Big\}\, \text{d}V

Developing the left-hand side first.

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q h_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(30)

= \left(\varepsilon \alpha_q \rho_q h_q u_q A\right)_{\text{e}} - \left(\varepsilon \alpha_q \rho_q h_q u_q A\right)_{\text{w}} - \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x} A\right)_{\text{e}} + \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x} A\right)_{\text{w}}

= F_{\text{e},q} h_{\text{e},q} + F_{\text{w},q} h_{\text{w},q} - \left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{e}} \left(h_{\text{E},q} - h_{\text{P},q}\right) + \left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{w}} \left(h_{\text{P},q} - h_{\text{W},q}\right)

Adopting upwind discretization for the enthalpy fluxes on the cell faces

$\begin{align} & F_{\text{e},q} h_{\text{e},q} = \text{max}\left(F_{\text{e},q},0\right) h_{\text{P},q} - \text{max}\left(-F_{\text{e},q},0\right) h_{\text{E},q}\\ & F_{\text{w},q} h_{\text{w},q} = \text{max}\left(F_{\text{w},q},0\right) h_{\text{W},q} - \text{max}\left(-F_{\text{w},q},0\right) h_{\text{P},q} \end{align}$

(31)

and rearranging terms

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q h_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(32)

= \left[ \text{max}\left(F_{\text{e},q},0\right) + \text{max}\left(-F_{\text{w},q},0\right) + \left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{e}} + \left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{w}} \right] h_{\text{P},q}

- \left[ \text{max}\left(-F_{\text{e},q},0\right) + \left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{e}} \right] h_{\text{E},q} - \left[ \text{max}\left(F_{\text{w},q},0\right) + \left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{w}} \right] h_{\text{W},q}

Denoting $\left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{e}} = D_{\text{e},q}$ and $\left(\varepsilon \alpha_q A \frac{\lambda_{\text{eff},q}}{c_{p,q}\delta x}\right)_{\text{w}} = D_{\text{w},q}$ and modifying the central coefficient slightly, since $\text{max}\left(F_{\text{e},q},0\right) = F_{\text{e},q} + \text{max}\left(-F_{\text{e},q},0\right)$ and $\text{max}\left(-F_{\text{w},q},0\right) = -F_{\text{w},q} + \text{max}\left(F_{\text{w},q},0\right)$ , Equation 32 reduces to

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q h_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(33)

= \left[ \text{max}\left(-F_{\text{e},q},0\right) + D_{\text{e},q} + \text{max}\left(F_{\text{w},q},0\right) + D_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q}

- \left[ \text{max}\left(-F_{\text{e},q},0\right) + D_{\text{e},q} \right] h_{\text{E},q} - \left[ \text{max}\left(F_{\text{w},q},0\right) + D_{\text{w},q} \right] h_{\text{W},q}

Identifying the coefficients of neighboring velocities as

$\begin{align} & a_{\text{e},q} = \text{max}\left(-F_{\text{e},q},0\right) + D_{\text{e},q}\\ & a_{\text{w},q} = \text{max}\left(F_{\text{w},q},0\right) + D_{\text{w},q} \end{align}$

(34)

leads to a general form

$\Leftrightarrow \iint\limits_{A}\!\!\!\!\!\!\!\!\!\!\!\subset\!\supset \left[\left(\varepsilon \alpha_q \rho_q h_q u_q\right)\cdot n - \left(\varepsilon \alpha_q \frac{\lambda_{\text{eff},q}}{c_{p,q}} \frac{\Delta h_q}{\delta x}\right)\cdot n\right]\, \text{d}A$

(35)

= \left[ a_{\text{e},q} + a_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}

Substituting the results of Equation 35 into Equation 29 and continuing with the development

$\Rightarrow \left[ a_{\text{e},q} + a_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}$

(36)

= \iiint\limits_V \Big\{ \left[\text{max}\left(\gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right]

-\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right) +\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right) -\nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right) + \nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) + s_{h,q} \Big\}\, \text{d}V

Integrating the right-hand side term by term. The phase change term reduces to

$\iiint\limits_V \left[\text{max}\left(\gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right] \text{d}V$

(37)

= \left[\text{max}\left(\Gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\Gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right]

The divergence terms, $-\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right)+\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right)$ , are handled equivalently to the divergence terms on the left-hand side of the equation

$\iiint\limits_V \left[ -\nabla\cdot\left(\tfrac{1}{2}\varepsilon \alpha_q \rho_q \left| u_q \right|^2 u_q\right)+\nabla\cdot\left(\varepsilon \alpha_q \rho_q \left( \mathbf{g} \cdot \mathbf{r} \right)u_q\right)\right] \text{d}V$

(38)

= -\tfrac{1}{2}\left( F_{\text{e},q} \left| u_q \right|^2_{\text{e}} - F_{\text{w},q} \left| u_q \right|^2_{\text{w}}\right) + \left[ F_{\text{e},q} \left( \mathbf{g} \cdot \mathbf{r} \right)_{\text{e}} - F_{\text{w},q} \left( \mathbf{g} \cdot \mathbf{r} \right)_{\text{w}}\right]

Heat source due to heat fluxes at the cell faces, $\nabla\cdot\left(\varepsilon \alpha_q \mathbf{q} \right)$ , are dropped out. Heat sources from an external source, such as neutronics solver, will be implemented later.

Friction work is replaced with the following simplified formulation

$\iiint\limits_V \left[\nabla\cdot\left(\varepsilon \alpha_q \boldsymbol{\tau} \cdot \mathbf{u}_q \right) \right] \text{d}V = \left( F_{\text{D}} + F_{\text{ID}} \right) u_{\text{P},q}$

(39)

where the drag force, $F_{\text{D}}$ , and the interfacial drag force, $F_{\text{ID}}$ , are calculated from the following relations

$F_{\text{D}} = -\tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q$

(40)

$F_{\text{ID}} = -K \left( u_q - u_p \right)$

(41)

Volume integral of volumetric enthalpy source to phase q is simply

$\iiint\limits_V s_{h,q} \text{d}V = S_{h,q}$

(42)

Combining equations Equation 36, Equation 37, Equation 38, Equation 39 and Equation 42 leads to

$\Rightarrow \left[ a_{\text{e},q} + a_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}$

(43)

= \left[\text{max}\left(\Gamma_{pq},0\right) \left(h_{\text{P,sat,}q} + \tfrac{1}{2} \left| u_p \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right) - \text{max}\left(-\Gamma_{pq},0\right) \left(h_{\text{P,}q} + \tfrac{1}{2} \left| u_q \right|^2_{\text{P}} - \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}}\right)\right]

-\tfrac{1}{2}\left( F_{\text{e},q} \left| u_q \right|^2_{\text{e}} - F_{\text{w},q} \left| u_q \right|^2_{\text{w}}\right) + \left[ F_{\text{e},q} \left( \mathbf{g} \cdot \mathbf{r} \right)_{\text{e}} - F_{\text{w},q} \left( \mathbf{g} \cdot \mathbf{r} \right)_{\text{w}}\right] + \left( F_{\text{D}} + F_{\text{ID}} \right) u_{\text{P},q} + S_{h,q}

Rearranging terms

$\Leftrightarrow \left[ a_{\text{e},q} + a_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}$

(44)

= -\tfrac{1}{2} \left[ F_{\text{e},q} \left| u_q \right|^2_{\text{e}} - F_{\text{w},q} \left| u_q \right|^2_{\text{w}} + \text{max}\left(-\Gamma_{pq},0\right) \left| u_q \right|^2_{\text{P}} - \text{max}\left(\Gamma_{pq},0\right) \left| u_p \right|^2_{\text{P}} \right]

+ \left[ F_{\text{e},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{e}} - F_{\text{w},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{w}} - \Gamma_{pq} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}} \right] + \left[ \text{max}\left(\Gamma_{pq},0\right) h_{\text{P,sat,}q} - \text{max}\left(-\Gamma_{pq},0\right) h_{\text{P},q}\right]

+ \left( F_{\text{D}} + F_{\text{ID}} \right) u_{\text{P},q} + S_{h,q}

At this point, the conservation of mass equation times cell specific enthalpy, $h_{\text{P},q}$ , is subtracted from Equation 44

$\Rightarrow \left[ a_{\text{e},q} + a_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) - \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}$

(45)

= -\tfrac{1}{2} \left[ F_{\text{e},q} \left| u_q \right|^2_{\text{e}} - F_{\text{w},q} \left| u_q \right|^2_{\text{w}} + \text{max}\left(-\Gamma_{pq},0\right) \left| u_q \right|^2_{\text{P}} - \text{max}\left(\Gamma_{pq},0\right) \left| u_p \right|^2_{\text{P}} \right]

+ \left[ F_{\text{e},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{e}} - F_{\text{w},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{w}} - \Gamma_{pq} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}} \right] + \left[ \text{max}\left(\Gamma_{pq},0\right) h_{\text{P,sat,}q} - \text{max}\left(-\Gamma_{pq},0\right) h_{\text{P},q}\right]

+ \left( F_{\text{D}} + F_{\text{ID}} \right) u_{\text{P},q} + S_{h,q} - \left( \Gamma_{pq} + S_{\text{mass},q} \right) h_{\text{P},q}

Splitting the mass source terms into negative and positive parts

$\Leftrightarrow \left[ a_{\text{e},q} + a_{\text{w},q} + \left( F_{\text{e},q} - F_{\text{w},q} \right) - \left( F_{\text{e},q} - F_{\text{w},q} \right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}$

(46)

= -\tfrac{1}{2} \left[ F_{\text{e},q} \left| u_q \right|^2_{\text{e}} - F_{\text{w},q} \left| u_q \right|^2_{\text{w}} + \text{max}\left(-\Gamma_{pq},0\right) \left| u_q \right|^2_{\text{P}} - \text{max}\left(\Gamma_{pq},0\right) \left| u_p \right|^2_{\text{P}} \right]

+ \left[ F_{\text{e},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{e}} - F_{\text{w},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{w}} - \Gamma_{pq} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}} \right] + \left[ \text{max}\left(\Gamma_{pq},0\right) h_{\text{P,sat,}q} - \text{max}\left(-\Gamma_{pq},0\right) h_{\text{P},q}\right]

+ \left( F_{\text{D}} + F_{\text{ID}} \right) u_{\text{P},q} + S_{h,q} - \left[ \text{max}\left(\Gamma_{pq},0\right) - \text{max}\left(-\Gamma_{pq},0\right) \right] h_{\text{P},q} - \left[ \text{max}\left(S_{\text{mass},q},0\right) - \text{max}\left(-S_{\text{mass},q},0\right)\right] h_{\text{P},q}

and rearranging the terms leads to the final form

$\Leftrightarrow \left[ a_{\text{e},q} + a_{\text{w},q} + \text{max}\left(\Gamma_{pq},0\right) + \text{max}\left(S_{\text{mass},q},0\right) + \text{max}\left(-\Gamma_{pq},0\right) \right] h_{\text{P},q} - a_{\text{e},q} h_{\text{E},q} - a_{\text{w},q} h_{\text{W},q}$

(47)

= \left[ \text{max}\left(-\Gamma_{pq},0\right) + \text{max}\left(-S_{\text{mass},q},0\right) \right] h_{\text{P},q} + \text{max}\left(\Gamma_{pq},0\right) h_{\text{P,sat},q}

- \tfrac{1}{2} \left[ F_{\text{e},q} \left| u_q \right|^2_{\text{e}} - F_{\text{w},q} \left| u_q \right|^2_{\text{w}} + \text{max}\left(-\Gamma_{pq},0\right) \left| u_q \right|^2_{\text{P}} - \text{max}\left(\Gamma_{pq},0\right) \left| u_p \right|^2_{\text{P}} \right]

+ \left[ F_{\text{e},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{e}} - F_{\text{w},q} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{w}} - \Gamma_{pq} \left(\mathbf{g} \cdot \mathbf{r}\right)_{\text{P}} \right] + \left( F_{\text{D}} + F_{\text{ID}} \right) u_{\text{P},q} + S_{h,q}

where:

$\alpha_q$		volume fraction of phase q [-],
$\gamma_{pq}$		volumetric mass transfer rate from phase p to q [kg/m³s],
$\Gamma_{pq}$		mass transfer rate from phase p to q [kg/s],
$\delta x$		distance between two adjacent cell centers [m],
$\varepsilon$		porosity, fluid fraction of cell volume [-],
$\lambda_{\text{eff},q}$		effective thermal conductivity of phase q [W/mK],
$\rho_q$		density of phase q [kg/m³],
$A$		face area [m²],
$C$		inertial loss coefficient [1/m],
$F$		face mass flow rate [kg/s],
$F_{\text{D}}$		drag force [N],
$F_{\text{ID}}$		interfacial drag force [N],
$\mathbf{g}$		acceleration of gravity (vector) [m/s²],
$h_q$		specific enthalpy of phase q [J/kg],
$h_{0,q}$		total enthalpy of phase q [J/kg],
$K$		interphase drag coefficient [kg/s],
$n$		face normal pointing out of the cell (1 for east face, -1 for west face),
$p$		pressure [Pa],
$\mathbf{r}$		location vector [m],
$s_{h,q}$		volumetric energy source to phase q [W/m³],
$S_{h,q}$		energy source to phase q [W],
$s_{\text{mass},q}$		volumetric mass source to phase q [kg/m³s],
$S_{\text{mass},q}$		mass source to phase q [kg/s],
$u_q$		velocity of phase q [m/s],
$V$		cell volume [m³],
subscript c1		refers to first cell,
subscript c2		refers to second cell,
subscript c3		refers to third cell,
subscript e		east face (the face in the negative x-direction),
subscript E		east neighbor (the cell in the negative x-direction),
subscript w		west face (the face in the positive x-direction),
subscript W		west neighbor (the cell in the positive x-direction),
subscript P		central cell (the cell for which the equation is written),
subscript q		current phase for which the equation is written (1 = primary, 2 = secondary),
subscript p		the other phase (p = 2 for q = 1 and p = 1 for q = 2),
subscript pq		indicates exchange between phases (e.g. $\Gamma_{pq}$ mass transfer from phase p to q).

Same as the conservation of momentum equations above, Equation 47 is solved in a phase-coupled manner, even though, due to the chosen interphase heat and mass transfer concept, the phase coupling coefficients are zero. As a reminder, the terms $\left[ \text{max}\left(\Gamma_{pq},0\right) + \text{max}\left(S_{\text{mass},q},0\right) \right] h_{\text{P},q}$ on the left-hand side and $\left[ \text{max}\left(-\Gamma_{pq},0\right) + \text{max}\left(-S_{\text{mass},q},0\right) \right] h_{\text{P},q}$ on the right-hand side result from the derivation procedure and are not considered as source terms. The terms $\text{max}\left(-\Gamma_{pq},0\right) h_{\text{P},q}$ on the left-hand side and $\text{max}\left(\Gamma_{pq},0\right) h_{\text{P,sat},q}$ on the right-hand side constitute the enthalpy source due to phase change. Note, however, that in the event of mass transfer from phase q, the terms $\text{max}\left(-\Gamma_{pq},0\right) h_{\text{P},q}$ on both sides of the equation cancel out.

Face Velocity Equation

The basic intent here is to find a relation for the face velocities, $u_{f,q}$ , that depend on the momentum solution of the adjacent cells and the pressure difference over the face in question. The importance of this formulation is paramount, for it determines whether a solution that satisfies continuity can coexist with a monotonous velocity field or not.

The first attempts at a solution algorithm based on a collocated formulation (both pressures and velocities are stored at cell centers) were found to produce highly oscillatory pressure and velocity fields. One of the first, or at least the most widely known, remedies, the Pressure Weighted Interpolation Method (PWIM), was proposed by Rhie and Chow ^[1] for single-phase flow. It quickly gained popularity in combination with the SIMPLE^[2] based methods. The original form of PWIM, however, is not without problems and several authors have since suggested modifications to it. The under-relaxation and time-step dependencies were first reported by Majumdar^[3] and Choi^[4], respectively. Recently, a well-written article was introduced by Pascau^[5], which includes a thorough and elegant derivation for the face velocity equation that deals with both the under-relaxation and time-step dependencies. In addition, the possible choices made in the interpolation of the terms in the face velocity equations and the different methods that result are well illustrated.^[5]

In the following, the derivation by Pascau^[5] is presented and slightly expanded for phase q. Although not directly evident, the resulting form is under-relaxation factor independent.

Starting from a compacted form of Equation 19, where all the other source terms, except the pressure gradient, are included in $S_{\text{mom,P},q}$ .

$a_{\text{P}} u_{\text{P},q} = \sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom,P},q} - \left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}}$

(48)

Under-relax

$\frac{1}{\alpha_u} a_{\text{P}} u_{\text{P},q} = \sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom,P},q} - \left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}} + \left( 1 - \alpha_u \right) \frac{1}{\alpha_u} a_{\text{P}} u_{\text{P},q}^*$

(49)

Multiply by $\alpha_u$

$a_{\text{P}} u_{\text{P},q} = \alpha_u \sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + \alpha_u S_{\text{mom,P},q} - \alpha_u \left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}} + \left( 1 - \alpha_u \right) a_{\text{P}} u_{\text{P},q}^*$

(50)

Divide by $a_{\text{P}}$

$u_{\text{P},q} = \alpha_u \frac{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom,P},q}} {a_{\text{P}}} - \alpha_u \frac{\left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}}} {a_{\text{P}}} + \left( 1 - \alpha_u \right) u_{\text{P},q}^*$

(51)

Writing a fictional momentum equation for face f that is identical in form to Equation 51

$u_{f,q} = \alpha_u \frac{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom},f,q}} {a_f} - \alpha_u \frac{\left( \varepsilon \alpha_q V \nabla p \right)_f} {a_f} + \left( 1 - \alpha_u \right) u_{f,q}^*$

(52)

and for momentum interpolated over face f

$\overline{u_{f,q}} = \alpha_u \frac{\overline{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom},f,q}}} {a_f} - \alpha_u \frac{\overline{\left( \varepsilon \alpha_q V \nabla p \right)_f}} {a_f} + \left( 1 - \alpha_u \right) \overline{u_{f,q}^*}$

(53)

Pressure Weighted Interpolation Method (PWIM):

$\frac{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom},f,q}} {a_f} = \frac{\overline{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom},f,q}}} {a_f}$

(54)

Solving $\frac{\overline{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom},f,q}}} {a_f}$ from Equation 53 and substituting it to Equation 52 leads to:

$u_{f,q} = \overline{u_{f,q}} + \alpha_u \left[ \frac{\overline{\left( \varepsilon \alpha_q V \nabla p \right)_f}} {a_f} - \frac{\left( \varepsilon \alpha_q V \nabla p \right)_f} {a_f} \right] + \left( 1 - \alpha_u \right) \left( u_{f,q}^* - \overline{u_{f,q}^*} \right)$

(55)

Phase-Coupled Face Velocity Equation

In this section, a phase-coupled version of the face velocity equations is developed. Similar ideas have been previously employed by Kunz et al. ^[6] and Vasquez and Ivanov ^[7]. Starting from a compacted form of Equation 19, where all the other source terms, except the pressure gradient, are included in $S_{\text{mom,P},q}$ .

$a_{\text{P}} u_{\text{P},q} = \sum_{nb=1}^N a_{\text{nb},q} u_{\text{nb},q} + S_{\text{mom,P},q} - \left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}}$

(56)

Under-relaxing the equation, as it has already been done in the solution of momentum equations in the previous step of the solution algorithm ( $\alpha_u$ = velocity under-relaxation factor). $u_{\text{P},q}^*$ denotes the value of $u_{\text{P},q}$ at the previous iteration.

$\Leftrightarrow \frac{1}{\alpha_u} a_{\text{P}} u_{\text{P},q} = \sum_{nb=1}^N a_{\text{nb},q} u_{\text{nb},q} + S_{\text{mom,P},q} - \left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}} + \left( 1 - \alpha_u \right) \frac{1}{\alpha_u} a_{\text{P}} u_{\text{P},q}^*$

(57)

Multiply by $\alpha_u$

$\Leftrightarrow a_{\text{P}} u_{\text{P},q} = \alpha_u \sum_{nb=1}^N a_{\text{nb},q} u_{\text{nb},q} + \alpha_u S_{\text{mom,P},q} - \alpha_u \left( \varepsilon \alpha_q V \nabla p \right)_{\text{P}} + \left( 1 - \alpha_u \right) a_{\text{P}} u_{\text{P},q}^*$

(58)

Casting Equation 58 into a vector form (over phases)

$\Rightarrow \begin{bmatrix} a_{\text{P},1} & 0 \\ 0 & a_{\text{P},2} \end{bmatrix} \begin{pmatrix} u_{\text{P},1}\\ u_{\text{P},2} \end{pmatrix} = \alpha_u \begin{pmatrix} \sum_{nb=1}^N a_{\text{nb},1} u_{\text{nb},1} + S_{\text{mom,P,1}}\\ \sum_{nb=1}^N a_{\text{nb},2} u_{\text{nb},2} + S_{\text{mom,P,2}} \end{pmatrix} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \begin{pmatrix} \alpha_{\text{P},1}\\ \alpha_{\text{P},2} \end{pmatrix} + \left( 1 - \alpha_u \right) \begin{bmatrix} a_{\text{P},1} & 0 \\ 0 & a_{\text{P},2} \end{bmatrix} \begin{pmatrix} u_{\text{P},1}^*\\ u_{\text{P},2}^* \end{pmatrix}$

(59)

Denoting $\mathbf{A} = \begin{bmatrix} a_{\text{P},1} & 0 \\ 0 & a_{\text{P},2} \end{bmatrix}$ , $\overrightarrow{u} = \begin{pmatrix} u_{\text{P},1} \\ u_{\text{P},2} \end{pmatrix}$ , $\overrightarrow{\alpha} = \begin{pmatrix} \alpha_{\text{P},1} \\ \alpha_{\text{P},2} \end{pmatrix}$ , and $\overrightarrow{u}^* = \begin{pmatrix} u_{\text{P},1}^* \\ u_{\text{P},2}^* \end{pmatrix}$ to simplify the notation.

$\Leftrightarrow \mathbf{A} \overrightarrow{u} = \alpha_u \begin{pmatrix} \sum_{nb=1}^N a_{\text{nb},1} u_{\text{nb},1} + S_{\text{mom,P,1}}\\ \sum_{nb=1}^N a_{\text{nb},2} u_{\text{nb},2} + S_{\text{mom,P,2}} \end{pmatrix} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \overrightarrow{\alpha} + \left( 1 - \alpha_u \right) \mathbf{A} \overrightarrow{u}^*$

(60)

Separating a part of the interphase friction from the source term

$\Leftrightarrow \mathbf{A} \overrightarrow{u} = \alpha_u \begin{pmatrix} \sum_{nb=1}^N a_{\text{nb},1} u_{\text{nb},1} + S_{\text{mom,P,1}}^*\\ \sum_{nb=1}^N a_{\text{nb},2} u_{\text{nb},2} + S_{\text{mom,P,2}}^* \end{pmatrix} + \alpha_u \begin{bmatrix} 0 & K \\ K & 0 \end{bmatrix} \overrightarrow{u} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \overrightarrow{\alpha} + \left( 1 - \alpha_u \right) \mathbf{A} \overrightarrow{u}^*$

(61)

and denoting $\overrightarrow{H} = \begin{pmatrix} \sum_{nb=1}^N a_{\text{nb},1} u_{\text{nb},1} + S_{\text{mom,P,1}}^*\\ \sum_{nb=1}^N a_{\text{nb},2} u_{\text{nb},2} + S_{\text{mom,P,2}}^* \end{pmatrix}$ , and $\mathbf{R} = \begin{bmatrix} 0 & K \\ K & 0 \end{bmatrix}$ .

$\Leftrightarrow \mathbf{A} \overrightarrow{u} = \alpha_u \overrightarrow{H} + \alpha_u \mathbf{R} \overrightarrow{u} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \overrightarrow{\alpha} + \left( 1 - \alpha_u \right) \mathbf{A} \overrightarrow{u}^*$

(62)

Subtracting $\mathbf{R} \overrightarrow{u}$ from both sides of the equation.

$\Leftrightarrow \left( \mathbf{A} - \mathbf{R} \right) \overrightarrow{u} = \alpha_u \overrightarrow{H} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \overrightarrow{\alpha} + \left( 1 - \alpha_u \right) \mathbf{A} \overrightarrow{u}^* - \left( 1 - \alpha_u \right) \mathbf{R} \overrightarrow{u}$

(63)

Since $\overrightarrow{u}$ is not available on the right-hand side, it is replaced by $\overrightarrow{u}^*$ .

$\Leftrightarrow \left( \mathbf{A} - \mathbf{R} \right) \overrightarrow{u} = \alpha_u \overrightarrow{H} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \overrightarrow{\alpha} + \left( 1 - \alpha_u \right) \left( \mathbf{A} - \mathbf{R} \right) \overrightarrow{u}^*$

(64)

Multiplying both sides by $\left( \mathbf{A} - \mathbf{R} \right)^{-1}$

$\Leftrightarrow \overrightarrow{u} = \alpha_u \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H} - \alpha_u \left( \varepsilon V \nabla p \right)_{\text{P}} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha} + \left( 1 - \alpha_u \right) \overrightarrow{u}^*$

(65)

Writing a fictitious momentum equation for face f that is identical in form to Equation 65

$\Rightarrow \overrightarrow{u_f} = \alpha_u \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f} - \alpha_u \left( \varepsilon V \nabla p \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f} + \left( 1 - \alpha_u \right) \overrightarrow{u_f}^*$

(66)

and for momentum interpolated over face f

$\Rightarrow \overline{\overrightarrow{u_f}} = \alpha_u \overline{\left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f}} - \alpha_u \overline{\left( \varepsilon V \nabla p \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}} + \left( 1 - \alpha_u \right) \overline{\overrightarrow{u_f}^*}$

(67)

When adopting the Pressure Weighted Interpolation Method (PWIM), the approach originally used by Rhie and Chow [1] can be interpreted as an interpolation of $\left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f}$ and $\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}$ over the face:

$\left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f} = \overline{\left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f}}, \qquad \left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f} = \overline{\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}}$

(68)

Solving $\overline{\left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f}}$ from Equation 67 and substituting it to Equation 66 leads to

$\Rightarrow \overrightarrow{u_f} = \overline{\overrightarrow{u_f}} + \alpha_u \left[\overline{\left( \varepsilon V \nabla p \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}} - \left( \varepsilon V \nabla p \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f} \right] + \left( 1 - \alpha_u \right) \left( \overrightarrow{u_f}^* - \overline{\overrightarrow{u_f}^*} \right)$

(69)

Rearranging variables and using $\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f} = \overline{\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}}$

$\Rightarrow \overrightarrow{u_f} = \overline{\overrightarrow{u_f}} + \alpha_u \left[\overline{\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f} \nabla p_f} - \overline{\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}} \nabla p_f \right] + \left( 1 - \alpha_u \right) \left( \overrightarrow{u_f}^* - \overline{\overrightarrow{u_f}^*} \right)$

(70)

Computation of face velocity involves the solution of two systems of equations for each face

$\begin{align} & \overrightarrow{x_{\text{c1}}} = \begin{pmatrix} x_{\text{c1},1} \\ x_{\text{c1},2} \end{pmatrix} = \begin{bmatrix} a_{\text{c1},1} & -K \\ -K & a_{\text{c1},2} \end{bmatrix}^{-1} \begin{pmatrix} \left( \varepsilon V \right)_{c1} \alpha_{\text{c1},1} \\ \left( \varepsilon V \right)_{c1} \alpha_{\text{c1},2} \end{pmatrix} \\ & \overrightarrow{x_{\text{c2}}} = \begin{pmatrix} x_{\text{c2},1} \\ x_{\text{c2},2} \end{pmatrix} = \begin{bmatrix} a_{\text{c2},1} & -K \\ -K & a_{\text{c2},2} \end{bmatrix}^{-1} \begin{pmatrix} \left( \varepsilon V \right)_{c2} \alpha_{\text{c2},1} \\ \left( \varepsilon V \right)_{c2} \alpha_{\text{c2},2} \end{pmatrix} \end{align}$

(71)

Linear interpolation of the first term inside the brackets, in Equation 70, yields

$\overline{\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f} \nabla p_f} = \beta \overrightarrow{x_{\text{c1}}} \nabla p_{\text{c1}} + \left( 1 - \beta \right) \overrightarrow{x_{\text{c2}}} \nabla p_{\text{c2}} = \overline{\overrightarrow{x_f} \nabla p_f}$

(72)

and the second term reduces to

$\overline{\left( \varepsilon V \right)_{f} \left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{\alpha_f}} \nabla p_f = \left[ \beta \overrightarrow{x_{\text{c1}}} + \left( 1 - \beta \right) \overrightarrow{x_{\text{c2}}} \right] \frac{p_{\text{c2}} - p_{\text{c1}}}{\delta x} = \frac{1}{\delta x} \overline{\overrightarrow{x_f}} \left( p_{\text{c2}} - p_{\text{c1}} \right)$

(73)

Substituting the results of Equation 72 and Equation 73 to Equation 70 leads to the final form, in which the equation is implemented in the code:

$\Rightarrow \overrightarrow{u_f} = \overline{\overrightarrow{u_f}} + \alpha_u \left[\overline{\overrightarrow{x_f} \nabla p_f} - \frac{1}{\delta x} \overline{\overrightarrow{x_f}} \left( p_{\text{c2}} - p_{\text{c1}} \right) \right] + \left( 1 - \alpha_u \right) \left( \overrightarrow{u_f}^* - \overline{\overrightarrow{u_f}^*} \right)$

(74)

Pressure & Volume Fraction Correction

The pressure & volume fraction correction method (pα-correction) starts with the idea of developing such corrections to pressure $p^'$ , phase volume fraction $\alpha_q^'$ , and in the process, to face mass flow rates $F_{f,q}^'$

$\begin{align} & p = p^* + p^' \\ & \alpha_q = \alpha_q^* + \alpha_q^', \qquad \qquad \sum{\alpha_q^'} = 0 \\ & F_{f,q} = F_{f,q}^* + F_{f,q}^' \end{align}$

(75)

that continuity of both phases (and consequently the sum of phases) is satisfied at the end of the iteration:

$\sum_{f=1}^N F_{f,q} = \left( \varepsilon \alpha_q \rho_q u_q A \right)_{\text{e}} - \left( \varepsilon \alpha_q \rho_q u_q A \right)_{\text{w}} = F_{\text{e},q} - F_{\text{w},q} = \Gamma_{pq} + S_{\text{mass},q}$

(76)

Starred values represent the latest available values and the dotted values stand for the corrections. The sum of volume fraction corrections is zero, which for a two-phase system means that $\alpha_1^' = - \alpha_2^'$ .

Face mass flow rates are functions of both velocity and volume fraction evaluated at the face. To obtain the influence of both pressure correction (through changes in face velocity) and volume fraction correction on face mass flow rate, face mass flow rate is derivated with respect to velocity and volume fraction.

$\begin{align} F_{f,q}^' &= D \left( F_{f,q} \right) = D \left( \varepsilon \alpha_q \rho_q u_q A \right)_f = \left( \varepsilon \rho_q A \right)_f \left[ D \left( \alpha_q u_q^' \right)_f + D \left( \alpha_q^' u_q \right)_f \right] \\ &= \left( \varepsilon \alpha_q \rho_q A \right)_f u_{f,q}^' + \left( \varepsilon \rho_q u_{q} A \right)_f \alpha_{f,q}^' \end{align}$

(77)

The relation between pressure and velocity corrections is obtained from equation Equation 74 by replacing the pressures of the adjacent cells by the corresponding pressure corrections and disregarding all the other terms. Note that the pressure difference, $p_{\text{c2}} - p_{\text{c1}}$ , has been reversed to produce positive pressure coefficients, $c_f$ .

$\begin{align} & \overrightarrow{u_f}^' = \alpha_u \frac{1}{\delta x} \overline{\overrightarrow{x_f}} \left( p_{\text{c1}} - p_{\text{c2}} \right) = \overrightarrow{c_f} \left( p_{\text{c1}} - p_{\text{c2}} \right) \qquad \text{or} \\ & u_{f,q}^' = \alpha_u \frac{1}{\delta x} \overline{x_{f,q}} \left( p_{\text{c1}} - p_{\text{c2}} \right) = c_{f,q} \left( p_{\text{c1}} - p_{\text{c2}} \right) \end{align}$

(78)

The development of the method begins with the phase continuity equations.

$\sum_{f=1}^N F_{f,q} = \Gamma_{pq} + S_{\text{mass},q}$

(79)

Introducing the corrected face mass flow rates from Equation 75 into Equation 79.

$\Rightarrow \sum_{f=1}^N F_{f,q}^' = \Gamma_{pq} + S_{\text{mass},q} - \sum_{f=1}^N F_{f,q}^*$

(80)

Writing the face mass flow corrections, $F_{f,q}^'$ , with the help of derivatives developed in Equation 77.

$\Rightarrow \left( \varepsilon \alpha_q \rho_q A \right)_{\text{e}} u_{\text{e},q}^' - \left( \varepsilon \alpha_q \rho_q A \right)_{\text{w}} u_{\text{w},q}^' + \left( \varepsilon \rho_q u_{q} A \right)_{\text{e}} \alpha_{\text{e},q}^' - \left( \varepsilon \rho_q u_{q} A \right)_{\text{w}} \alpha_{\text{w},q}^' = \Gamma_{pq} + S_{\text{mass},q} - \sum_{f=1}^N F_{f,q}^*$

(81)

Substituting Equation 78 into Equation 81 to express the velocity correction in terms of adjacent pressure corrections.

$\Rightarrow \left( \varepsilon \alpha_q \rho_q A \right)_{\text{e}} c_{e,q} \left( p_{\text{P}}^' - p_{\text{E}}^' \right) - \left( \varepsilon \alpha_q \rho_q A \right)_{\text{w}} c_{w,q} \left( p_{\text{W}}^' - p_{\text{P}}^' \right)$

+ \left( \varepsilon \rho_q u_{q} A \right)_{\text{e}} \alpha_{\text{e},q}^' - \left( \varepsilon \rho_q u_{q} A \right)_{\text{w}} \alpha_{\text{w},q}^' = \Gamma_{pq} + S_{\text{mass},q} - \sum_{f=1}^N F_{f,q}^*

(82)

Denoting $a_{\text{e},p,q} = \left( \varepsilon \alpha_q \rho_q A \right)_{\text{e}} c_{e,q}$ , $a_{\text{w},p,q} = \left( \varepsilon \alpha_q \rho_q A \right)_{\text{w}} c_{w,q}$ , and $a_{\text{P},p,q} = a_{\text{e},p,q} + a_{\text{w},p,q}$ to simplify the notation:

$\Rightarrow a_{\text{P},p,q} p_{\text{P}}^' - a_{\text{e},p,q} p_{\text{E}}^' - a_{\text{w},p,q} p_{\text{W}}^'$

+ \left( \varepsilon \rho_q u_{q} A \right)_{\text{e}} \alpha_{\text{e},q}^' - \left( \varepsilon \rho_q u_{q} A \right)_{\text{w}} \alpha_{\text{w},q}^' = \Gamma_{pq} + S_{\text{mass},q} - \sum_{f=1}^N F_{f,q}^*

(83)

Upwind discretization is adopted for the remaining face values on the left-hand side:

$\begin{align} & \left( \varepsilon \rho_q u_{q} A \right)_{\text{e}} \alpha_{\text{e},q}^' = \left( \varepsilon \rho_q \right)_{\text{P}} A_{\text{e}} \text{max}\left( u_{\text{e},q},0 \right) \alpha_{\text{P},q}^' - \left( \varepsilon \rho_q \right)_{\text{E}} A_{\text{e}} \text{max}\left( -u_{\text{e},q},0 \right) \alpha_{\text{E},q}^' \\ & \left( \varepsilon \rho_q u_{q} A \right)_{\text{w}} \alpha_{\text{w},q}^' = \left( \varepsilon \rho_q \right)_{\text{W}} A_{\text{w}} \text{max}\left( u_{\text{w},q},0 \right) \alpha_{\text{W},q}^' - \left( \varepsilon \rho_q \right)_{\text{P}} A_{\text{w}} \text{max}\left( -u_{\text{w},q},0 \right) \alpha_{\text{P},q}^' \end{align}$

(84)

Denoting $a_{\text{e},\alpha,q} = \left( \varepsilon \rho_q \right)_{\text{E}} A_{\text{e}} \text{max}\left( -u_{\text{e},q},0 \right)$ , $a_{\text{w},\alpha,q} = \left( \varepsilon \rho_q \right)_{\text{W}} A_{\text{w}} \text{max}\left( u_{\text{w},q},0 \right)$ , and $a_{\text{P},\alpha,q} = \left( \varepsilon \rho_q \right)_{\text{P}} A_{\text{e}} \text{max}\left( u_{\text{e},q},0 \right) + \left( \varepsilon \rho_q \right)_{\text{P}} A_{\text{w}} \text{max}\left( -u_{\text{w},q},0 \right)$ and substituting the upwinded terms to Equation 83

$\Rightarrow a_{\text{P},p,q} p_{\text{P}}^' - a_{\text{e},p,q} p_{\text{E}}^' - a_{\text{w},p,q} p_{\text{W}}^'$

+ a_{\text{P},\alpha,q} \alpha_{\text{P},q}^' - a_{\text{e},\alpha,q} \alpha_{\text{E},q}^' - a_{\text{w},\alpha,q} \alpha_{\text{W},q}^' = \Gamma_{pq} + S_{\text{mass},q} - \sum_{f=1}^N F_{f,q}^*

(85)

Writing Equation 85 for both phases and substituting $\alpha_1^' = -\alpha_2^'$ to the primary phase equation leads to the final form:

$\begin{cases} a_{\text{P},p,1} p_{\text{P}}^' - a_{\text{e},p,1} p_{\text{E}}^' - a_{\text{w},p,1} p_{\text{W}}^' - a_{\text{P},\alpha,1} \alpha_{\text{P},2}^' + a_{\text{e},\alpha,1} \alpha_{\text{E},2}^' + a_{\text{w},\alpha,1} \alpha_{\text{W},2}^' = \Gamma_{21} + S_{\text{mass},1} - \sum_{f=1}^N F_{f,1}^* \\ a_{\text{P},p,2} p_{\text{P}}^' - a_{\text{e},p,2} p_{\text{E}}^' - a_{\text{w},p,2} p_{\text{W}}^' + a_{\text{P},\alpha,2} \alpha_{\text{P},2}^' - a_{\text{e},\alpha,2} \alpha_{\text{E},2}^' - a_{\text{w},\alpha,2} \alpha_{\text{W},2}^' = \Gamma_{12} + S_{\text{mass},2} - \sum_{f=1}^N F_{f,2}^* \end{cases}$

(86)

A single system of equations is formed from Equation 86 by stacking up the primary phase equations on top of the secondary phase equations. Pressure correction is chosen as the unknown for the primary phase equations and secondary phase volume fraction correction as the unknown for the secondary phase equations. The resulting system of equations (for a 1D channel consisting of three adjacent cells) is shown below.

$\begin{bmatrix} a_{\text{P,c1},p,1} & -a_{\text{e,c1},p,1} & 0 & -a_{\text{P,c1},\alpha,1} & a_{\text{e,c1},\alpha,1} & 0 \\ -a_{\text{w,c2},p,1} & a_{\text{P,c2},p,1} & -a_{\text{e,c2},p,1} & a_{\text{w,c2},\alpha,1} & -a_{\text{P,c2},\alpha,1} & a_{\text{e,c2},\alpha,1} \\ 0 & -a_{\text{w,c3},p,1} & a_{\text{P,c3},p,1} & 0 & a_{\text{w,c3},\alpha,1} & -a_{\text{P,c3},\alpha,1} \\ a_{\text{P,c1},p,2} & -a_{\text{e,c1},p,2} & 0 & a_{\text{P,c1},\alpha,2} & -a_{\text{e,c1},\alpha,2} & 0 \\ -a_{\text{w,c2},p,2} & a_{\text{P,c2},p,2} & -a_{\text{e,c2},p,2} & -a_{\text{w,c2},\alpha,2} & a_{\text{P,c2},\alpha,2} & -a_{\text{e,c2},\alpha,2} \\ 0 & -a_{\text{w,c3},p,2} & a_{\text{P,c3},p,2} & 0 & -a_{\text{w,c3},\alpha,2} & a_{\text{P,c3},\alpha,2} \end{bmatrix} \begin{pmatrix} p_{\text{c1}}^' \\ p_{\text{c2}}^' \\ p_{\text{c3}}^' \\ \alpha_{\text{c1,2}}^'\\ \alpha_{\text{c2,2}}^'\\ \alpha_{\text{c3,2}}^' \end{pmatrix} = \begin{pmatrix} b_{\text{c1,1}}\\ b_{\text{c2,1}}\\ b_{\text{c3,1}}\\ b_{\text{c1,2}}\\ b_{\text{c2,2}}\\ b_{\text{c3,2}} \end{pmatrix}$

(87)

For some linear solvers it might be beneficial to load the equations in a different order. The ordering shown in Equation 21 for the coupled momentum equations produces a significant reduction in bandwidth for large systems of equations, but for such a small number of cells (~20 cells in a single channel) the efficiency gain is probably insignificant.

After the solution of Equation 87 face velocities are corrected with

$u_{f,q} = u_{f,q}^* + \alpha_u \frac{1}{\delta x} \overline{x_{f,q}} \left( p_{\text{c1}} - p_{\text{c2}} \right)$

(88)

phase volume fractions are corrected with

$\begin{align} & \alpha_1 = \alpha_1^* - \alpha_2^' \\ & \alpha_2 = \alpha_2^* + \alpha_2^' \end{align}$

(89)

and pressure is corrected with

$p = p^* + \alpha_p p^'$

(90)

Note that the pressure correction can be under-relaxed, with an under-relaxation factor $\alpha_p$ , while the full amount of the correction needs to be applied to the face velocities and volume fractions in order to produce face mass flow rates that satisfy continuity. Cell-centered velocities can still be under-relaxed in the usual manner when solving the momentum equations, since under-relaxation has been taken into account in the face velocity equations and the pressure & volume fraction correction equations.

Solution Algorithm

File:SolutionAlgorithmSmall.png

Figure 4: Solution algorithm

The numerical solution of the governing equations relies on a pressure-based segregated approach, which means that pressure has been selected as one of the solution variables, instead of density, and each solution variable is solved from its own system of equations, one after the other. The only exception, which moves the classification of the solution method towards coupled methods, is that pressure and volume fraction corrections are solved together from a single system of equations.

The implemented solution algorithm belongs to the SIMPLE family of algorithms, extended to multiphase flow from the original single-phase SIMPLE algorithm of Patankar^[2], and bears close resemblance to both the Coupled Phasic Exchange algorithm of Kunz et al. ^[6] and the Phase Coupled SIMPLE algorithm of Vasquez and Ivanov ^[7]. A schematic of the solution algorithm is shown in Figure 4.

Momentum Transfer

The forces that drive the fluid flow are placed on the right-hand side of the momentum conservation equations, best seen in Equation 16 before the linearization and other modifications. These are called the momentum source terms and they include the effects caused by phase change, pressure gradient, gravitational forces, pressure loss caused by the porous material, interfacial friction, and possibly user-defined momentum sources:

$S = \left[\text{max}\left(\Gamma_{pq},0\right)u_p - \text{max}\left(-\Gamma_{pq},0\right)u_q\right] - \varepsilon \alpha_q V \nabla p + \varepsilon \alpha_q \rho_q V g - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q - K\left(u_q - u_p\right) + S_{\text{mom},q}$

(91)

At this point a parameter fundamental to internal flows is introduced, since it often appears in 1D pressure drop correlations. Hydraulic diameter $D_{\text{h}}$ [m] is defined in terms of either the flow area, $A_{\text{f}}$ [m²], and wetted perimeter, $P_{\text{w}}$ [m], or the fluid volume, $V_{\text{f}}$ [m³], and wetted surface area, $A_{\text{w}}$ [m²], as follows

$D_{\text{h}} = \frac{4 A_{\text{f}}}{P_{\text{w}}} = \frac{4 V_{\text{f}}}{A_{\text{w}}} = \frac{4 \varepsilon V }{A_{\text{w}}}$

(92)

where fluid volume, $V_{\text{f}}$ [m³], can be expressed as the product of total volume, $V$ [m³], and porosity, $\varepsilon$ [-].

Pressure loss

The solid structures inside the flow channel exert a force that resists the fluid flow. The pressure loss (or drag) is modelled as an inertial resistance:

$F_{\text{D}} = - \tfrac{1}{2} \varepsilon^3 \alpha_q \rho_q C V \left| u_q \right| u_q$

(93)

The inertial loss coefficient, $C$ [1/m], is defined as the ratio of the moody friction factor to hydraulic diameter:

$C = \frac{f}{D_{\text{h}}}$

(94)

The friction factor, $f$ [-], depends on the phasic Reynolds number:

${\text{Re}}_{q} = \frac{\rho_q |u_{q}| D_{\text{t}}}{\mu_{q}}$

(95)

If the flow is laminar, the moody friction factor for laminar flow is used, and if the flow is turbulent, the friction factor is calculated according to the correlation of Petukhov.

$f = \begin{cases} 64 \, / \, {\text{Re}}_{q} \\ \\ \left[ 0.79 \, \text{ln} \left({\text{Re}}_{q}\right) - 1.64 \right]^{-2} \end{cases} \begin{align} & {\text{Re}}_{q} < 2300 \\ & \\ & {\text{Re}}_{q} \ge 2300 \end{align}$

(96)

Interface momentum transfer

$F_{\text{ID}} = - K\left(u_q - u_p\right)$

(97)

Heat and Mass Transfer

The chosen heat transfer model is a combination of a well-known wall boiling model (RPI^[8]) and an interfacial heat transfer model to account for the mass transfer between the phases. Some minor additions have been made in the implementation of the wall boiling model to extend applicability towards dry-out conditions. However, when approaching dry-out conditions, the results can not be considered accurate, since the model assumptions become unrealistic and the sub-model correlations go out of their valid range. In such occasions the results should be compared against a critical heat flux correlation for the specified geometry.

To facilitate coupling with fuel behavior codes, two types of thermal boundary conditions were created: either a specified heat flux or a surface temperature distribution can be provided, and the other will be calculated.

To make a distinction between the hydraulic diameter, that includes all the surfaces of the flow channel and is used to calculate the pressure drop along the length of the channel, a new parameter, called thermal diameter $D_{\text{t}}$ [m], is defined analogously to the hydraulic diameter above. Instead of using the wetted perimeter, $P_{\text{w}}$ [m], and wetted surface area, $A_{\text{w}}$ [m²], the thermal diameter is defined using the heated perimeter, $P_{\text{h}}$ [m], and heated surface area, $A_{\text{h}}$ [m²]:

$D_{\text{t}} = \frac{4 A_{\text{f}}}{P_{\text{h}}} = \frac{4 V_{\text{f}}}{A_{\text{h}}} = \frac{4 \varepsilon V }{A_{\text{h}}}$

(98)

Surface heat transfer

The RPI wall boiling model is a mechanistic model developed at Rensselaer Polytechnic Institute ^[8], in which the wall heat flux has been divided into three separate parts: single-phase convective heat flux into liquid outside the influence of the nucleation sites (convection), convective heat flux enhanced by the detaching bubbles within the influence of the nucleation sites (quenching) and evaporation heat transfer directly into vapor through a superheated liquid microlayer (evaporation). The main idea of mechanistic boiling models is to identify the separate mechanisms that affect the heat transfer and to handle them separately, in order to better account for the various multidimensional phenomena related to two-phase flows. Such models should therefore be less dependent, for example, on flow conditions and the geometry being modeled.

The RPI boiling model was originally developed for the subcooled boiling regime inside vertical tubes. With certain modifications, such as the choice of model correlations and the addition of heat flux partitioning, it has been extended to nucleate boiling regime, and even conditions approaching dry-out.

In the original model [8] the wall heat flux, $q_{\text{w}}$ [W/m²], is partitioned into three parts, as follows

$q_{\text{w}} = q_{\text{l}} + q_{\text{q}} + q_{\text{e}}$

(99)

where:

$q_{\text{l}}$		convective heat flux to the liquid phase outside the zone of influence of the detaching bubbles [W/m²],
$q_{\text{q}}$		convective heat flux to the liquid phase enhanced by detaching bubbles (quenching) [W/m²] and
$q_{\text{e}}$		evaporation heat flux directly to the vapor phase [W/m²].

The wall surface area is divided into two fractions: an area fraction affected, $A_2$ [-], and unaffected, $A_1$ [-], by the departure of the bubbles.

$A_1 + A_2 = 1$

(100)

The area fraction influenced by boiling, $A_2$ [-], is calculated from bubble departure diameter, $d_{\text{dep}}$ [m], nucleation site density, $N$ [1/m²], and a correction proposed by Del Valle and Kenning ^[9], $K$ [-]. It is limited to less than unity.

$\begin{align} & A_2 = \text{min} \left( 1, K \frac{\pi d_{\text{dep}}^2 N}{4} \right) \\ & K = 4.8 e^{-Ja / 80} \\ & Ja = \frac{\rho_1 c_{p,1} \left( T_{\text{sat}} - T_1 \right)}{\rho_2 \left( h_{2,\text{sat}} - h_1 \right)} \end{align}$

(101)

The area fraction unaffected by boiling is then calculated by subtracting the affected area fraction from unity, but due to numerical reasons a lower limit of 10^-4 has been set.

$A_1 = \text{max} \left( 10^{-4}, 1 - A_2 \right)$

(102)

The convective heat flux to the liquid phase can then be calculated from the area fraction unaffected by boiling, $A_1$ [-], surface temperature, $T_{\text{w}}$ [K], bulk liquid temperature, $T_{1}$ [K], and the convective heat transfer coefficient, $h_{\text{c},1}$ [W/m²K].

$q_{\text{l}} = A_1 h_{\text{c},1} \left( T_{\text{w}} - T_{1}\right)$

(103)

Here the convective heat transfer coefficient, $h_{\text{c},1}$ [W/m²K], is calculated based on the Dittus and Boelter equation, as introduced by McAdams^[10].

$h_{\text{c},1} = 0.023 \text{Re}_1^{0.8} \text{Pr}_1^{0.4} \frac{\lambda_1}{D_{\text{t}}}$

(104)

The phasic Reynolds and Prandtl numbers are defined, respectively, as follows

${\text{Re}}_{q} = \frac{\rho_q |u_{q}| D_{\text{t}}}{\mu_{q}}$

(105)

${\text{Pr}}_{q} = \frac{\mu_{q} c_{p,q}}{\lambda_{q}}$

(106)

where $\lambda_{q}$ thermal conductivity of phase q [W/m K], $\mu_{q}$ dynamic viscosity of phase q [Ns/m²], $\rho_q$ density of phase q [kg/m³], $c_{p,q}$ specific isobaric heat capacity of phase q [J/kgK], $D_{\text{t}}$ thermal diameter [m], and $u_q$ velocity of phase q [m/s].

The quenching heat flux is calculated from an equation identical in form to Equation 103

$q_{\text{q}} = A_2 h_{\text{q}} \left( T_{\text{w}} - T_{1}\right)$

(107)

The quenching heat transfer coefficient, $h_{\text{q}}$ , is taken according to Del Valle and Kenning ^[9]:

$h_{\text{q}} = 2 \lambda_1 f_{\text{dep}} \sqrt{\frac{c_{p,1} \rho_1}{f_{\text{dep}} \pi \lambda_1}}$

(108)

The evaporation heat flux is the product of evaporation mass flux and the difference between saturated vapor enthalpy and liquid enthalpy:

$q_{\text{e}} = \dot{m} \left( h_{2,\text{sat}} - h_{1} \right)$

(109)

Here the evaporation mass flux, $\dot{m}$ [kg\m²s], can be written in terms of nucleation site density, $N$ [1\m²], bubble departure frequency, $f_{\text{dep}}$ [1/s], bubble departure diameter, $d_{\text{dep}}$ [m], and vapor density, $\rho_{2,\text{sat}}$ [kg\m³].

$\dot{m} = \frac{\pi d_{\text{dep}}^3}{6} \rho_{2,\text{sat}} f_{\text{dep}} N$

(110)

Next, the model correlations implemented for the terms appearing in the preceding equations are listed. Bubble departure frequency is calculated according to the presentation of Cole^[11]:

$f_{\text{dep}} = \sqrt{\frac{4 g \left( \rho_1 - \rho_{2,\text{sat}}\right)}{3 d_{\text{dep}} \rho_1}}$

(111)

Bubble departure diameter is taken according to Tolubinski and Kostanchuk ^[12]:

$d_{\text{dep}} = \text{max} \left[ \text{min} \left( d_{\text{ref}} \; e^{-\frac{T_{\text{sat}} - T_1}{{\Delta T}_{\text{ref}}}}, d_{\text{max}} \right), d_{\text{min}} \right]$

(112)

where ${\Delta T}_{\text{ref}}$ = 45 K, $d_{\text{ref}}$ = 0.6e^-3 m, $d_{\text{min}}$ = 1.0e^-6 m, and $d_{\text{max}}$ = 1.4e^-3 m.

Nucleation site density, as presented by Lemmert and Chawla ^[13]

$N = \left[ 210 \left( T_{\text{w}} - T_{\text{sat}} \right) \right]^{1.805}$

(113)

is slightly modified to remove the fractional powers of temperature

$\Rightarrow N = N_{\text{ref}} \left( \frac{T_{\text{w}} - T_{\text{sat}}} { {\Delta T}_{\text{ref}} } \right)^p$

(114)

where $N_{\text{ref}} = \left( 210 {\Delta T}_{\text{ref}} \right)^p$ , ${\Delta T}_{\text{ref}}$ = 10 K, and $p$ = 1.805.

Finally, an additional level of partitioning is introduced to extend applicability towards high void fractions:

$q_{\text{w}} = f_{\text{liquid}} \left( q_{\text{l}} + q_{\text{q}} + q_{\text{e}} \right) + \left( 1 - f_{\text{liquid}} \right) q_{\text{v}}$

(115)

where $f_{\text{liquid}}$ is a coefficient that defines the fraction of the heat flux that is transferred to the liquid phase (liquid fraction) [-] and $q_{\text{v}}$ is the convective heat flux to the vapor phase [W/m²].

The liquid fraction, $f_{\text{liquid}}$ , is calculated according to the presentation of Lavieville et al. ^[14]:

$\begin{array}{ll} \begin{cases} f_{\text{liquid}} = 1 - \frac{1}{2} e^{-20 \left( \alpha_1 - \alpha_{\text{crit}} \right)}, \\ f_{\text{liquid}} = \frac{1}{2} \left( \frac{\alpha_1}{\alpha_{\text{crit}}} \right)^{20 \, \alpha_{\text{crit}} }, \end{cases} & \begin{align} & \alpha_1 \ge \alpha_{\text{crit}} \\ & \alpha_1 < \alpha_{\text{crit}} \end{align} \end{array}$

(116)

where $\alpha_1$ is the liquid volume fraction [-] and $\alpha_{\text{crit}}$ is the critical liquid volume fraction provided by the user. Lavieville et al. proposed a critical liquid volume fraction of $\alpha_{\text{crit}} = 0.2$ .

Convective heat flux to the vapor phase is calculated similar to Equation 103

$q_{\text{v}} = h_{\text{c},2} \left( T_{\text{w}} - T_{2}\right)$

(117)

and the convective heat transfer coefficient, $h_{\text{c},2}$ [W/m²K], is also calculated based on the Dittus and Boelter equation [10].

$h_{\text{c},2} = 0.023 \text{Re}_2^{0.8} \text{Pr}_2^{0.4} \frac{\lambda_2}{D_{\text{t}}}$

(118)

The wall heat flux, in Equation 115, needs to be transformed into a form in which it can be casted as a source term to the relevant enthalpy equations. This is achieved by multiplying the surface heat fluxes by the ratio of heated surface area and volume, $A_{\text{h}} / V$ [1/m]. There are several ways the heat flux can be distributed between the phases (and the phase interface), but the chosen one is capable of producing vapor in subcooled conditions (subcooled boiling).

$\begin{align} & Q_{\text{l}} = f_{\text{liquid}} \left( q_{\text{l}} + q_{\text{q}} \right) \frac{A_{\text{h}}}{V} \\ & Q_{\text{e}} = f_{\text{liquid}} \, q_{\text{e}} \frac{A_{\text{h}}}{V} \\ & Q_{\text{v}} = \left( 1 - f_{\text{liquid}} \right) q_{\text{v}} \frac{A_{\text{h}}}{V} \end{align}$

(119)

Here, $Q_{\text{l}}$ [W/m³] is added to the liquid enthalpy equations, $Q_{\text{v}}$ [W/m³] is added to the vapor enthalpy equations, and $Q_{\text{e}}$ [W/m³] is added to the phase interface, as shown in the following subsection. Using the definition of the thermal diameter in Equation 98, the volumetric heat fluxes in Equation 119 can be written in a more useful form:

$\begin{align} & Q_{\text{l}} = f_{\text{liquid}} \left( q_{\text{l}} + q_{\text{q}} \right) \frac{4 \varepsilon}{D_{\text{t}}} \\ & Q_{\text{e}} = f_{\text{liquid}} \, q_{\text{e}} \frac{4 \varepsilon}{D_{\text{t}}} \\ & Q_{\text{v}} = \left( 1 - f_{\text{liquid}} \right) q_{\text{v}} \frac{4 \varepsilon}{D_{\text{t}}} \end{align}$

(120)

Interfacial heat transfer

The heat transfer coefficient from phase q to the interface between the two phases, $h_{q \text{i}}$ [W/m²K], is calculated according to the correlation of Ranz and Marshall ^[15]:

$h_{q \text{i}} = \left( 2 + 0.6 \text{Re}_{\text{r},q}^{\frac{1}{2}} \text{Pr}_{q}^{\frac{1}{3}} \right) \frac{\lambda_q}{d_{\text{b}}}$

(121)

where ${\text{Re}}_{\text{r},q}$ is the relative Reynolds number for phase q

${\text{Re}}_{\text{r},q} = \frac{\rho_q |u_{\text{r}}| d_{\text{b}}}{\mu_{q}}$

(122)

and $\text{Pr}_{q}$ is the phasic Prandtl number in Equation 106, $\lambda_q$ is the thermal conductivity of phase q [W/m K], $\rho_q$ is the density of phase q [kg/m³], $\mu_{q}$ dynamic viscosity of phase q [Ns/m²], $d_{\text{b}}$ is the bubble diameter [m], and $u_{\text{r}}$ is the velocity difference between the phases, $u_2 - u_1$ [m/s].

The interfacial volumetric heat fluxes from phase q to the interface between the phases, $Q_{q \text{i}}$ [W/m³], are calculated from the following relation:

$Q_{q \text{i}} = h_{q \text{i}} \left( T_q - T_{\text{sat}} \right) \frac{A_{\text{i}}}{V} - \frac{{\text{d}m}_q}{\text{d}t} h_{\text{int},q}$

(123)

where $h_{q \text{i}}$ is the interfacial heat transfer coefficient [W/m²K], $T_q$ is the temperature of phase q [K], $T_{\text{sat}}$ is the saturation temperature [K], $A_{\text{i}} / V$ is the interface area density [1/m], ${\text{d}m}_q / \text{d}t$ is the mass transfer rate to phase q [kg/m³s], and $h_{\text{int},q}$ is the enthalpy of phase q at the interface [J/kg]. The last term on the right-hand side is due to phase change.

A rough approximation for the interface area density is obtained by assuming spherical bubbles and droplets with a constant diameter, $d_{\text{b}}$ [m], and assuming that the flow regime changes from bubbly flow to droplet flow at a vapor fraction of $\alpha_1 = \alpha_2 = 0.5$ .

$\frac{A_{\text{i}}}{V} = \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}}$

(124)

When Equation 124 is combined with Equation 123, the interfacial volumetric heat fluxes are reduced to:

$Q_{q \text{i}} = h_{q \text{i}} \left( T_q - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} - \frac{{\text{d}m}_q}{\text{d}t} h_{\text{int},q}$

(125)

Mass transfer

By enforcing energy balance on the phase interface, i.e. setting the sum of heat fluxes to the interface equal to zero, allows for the mass transfer rate to be solved from the resulting equation.

$Q_{1 \text{i}} + Q_{2 \text{i}} + Q_{\text{e}} = 0$

(126)

$\Leftrightarrow h_{1 \text{i}} \left( T_1 - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} - \frac{{\text{d}m}_1}{\text{d}t} h_{\text{int},1} + h_{2 \text{i}} \left( T_2 - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} - \frac{{\text{d}m}_2}{\text{d}t} h_{\text{int},2} + f_{\text{liquid}} \, q_{\text{e}} \frac{4 \varepsilon}{D_{\text{t}}} = 0$

(127)

Denoting volumetric evaporation rate [kg/m³s], i.e. volumetric mass transfer rate from phase 1 to phase 2, as $\gamma_{12} = {\text{d}m}_2 / \text{d}t = - {\text{d}m}_1 / \text{d}t$ leads to

$\Leftrightarrow h_{1 \text{i}} \left( T_1 - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} + h_{2 \text{i}} \left( T_2 - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} + f_{\text{liquid}} \, q_{\text{e}} \frac{4 \varepsilon}{D_{\text{t}}} = \gamma_{12} h_{\text{int},2} - \gamma_{12} h_{\text{int},1}$

(128)

$\Leftrightarrow \gamma_{12} = \frac{ h_{1 \text{i}} \left( T_1 - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} + h_{2 \text{i}} \left( T_2 - T_{\text{sat}} \right) \frac{6 \text{min} \left( \alpha_1, \alpha_2 \right)}{d_{\text{b}}} + f_{\text{liquid}} \, q_{\text{e}} \frac{4 \varepsilon}{D_{\text{t}}} } {h_{\text{int},2} - h_{\text{int},1}}$

(129)

The interfacial enthalpies are upwinded based on the direction of mass transfer:

$\begin{array}{lll} \begin{cases} h_{\text{int},1} = h_1, \\ h_{\text{int},1} = h_{\text{sat},1}, \end{cases} & \begin{align} & h_{\text{int},2} = h_{\text{sat},2}, \\ & h_{\text{int},2} = h_2, \end{align} & \begin{align} & \gamma_{12} \ge 0 \\ & \gamma_{12} < 0 \end{align} \end{array}$

(130)

References

↑ Rhie, C. M., Chow, W. L., Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal 21 (11), pp. 1525-1532, 1983.
↑ ^2.0 ^2.1 Patankar, S. V., Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation. 1980.
↑ Majumdar, S., Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids. Numerical Heat Transfer 13 (1), pp. 125-132, 1988.
↑ Choi, S. K., Note on the use of momentum interpolation method for unsteady flows. Numerical Heat Transfer: Part A: Applications 36 (5), pp. 545-550, 1999.
↑ ^5.0 ^5.1 ^5.2 Pascau, A., Cell face velocity alternatives in a structured colocated grid for the unsteady navier-stokes equations. International Journal for Numerical Methods in Fluids 65 (7), pp. 812-833, 2011.
↑ ^6.0 ^6.1 Kunz, R. F., Siebert, B. W., Cope, W. K., Foster, N. F., Antal, S. P., Ettorre, S. M., A coupled phasic exchange algorithm for three–dimensional multi–field analysis of heated flows with mass transfer. Computers & Fluids 27 (7), pp. 741–768, 1998.
↑ ^7.0 ^7.1 Vasquez, S. A., Ivanov, V. A., A phase coupled method for solving multiphase problems in unstructured meshes. In: Proceedings of ASME FEDSM’00: ASME 2000 fluids engineering division summer meeting, Boston. pp. 743-748, 2000.
↑ ^8.0 ^8.1 Kurul, N., Podowski, M. Z., On the Modeling of Multidimensional Effects in Boiling Channels. In: ANS Proceedings of 1991 National Heat Transfer Conference, Minneapolis, Minnesota, Vol. 5. pp. 30-40, 1991. ISBN: 0-89448-162-2
↑ ^9.0 ^9.1 Del Valle M., V. H., Kenning, D. B. R., Subcooled flow boiling at high heat flux. Int. J. Heat Mass Transfer, Vol. 28, No. 10, pp.1907-1920, 1985.
↑ McAdams, W. H., Heat Transmission, 2nd Edition, McGraw-Hill, New York, 1942.
↑ Cole, R., A photographic study of pool boiling in the region of CHF. AIChEJ, 6 pp. 533-542, 1960.
↑ Tolubinski, V. I., Kostanchuk, D. M., Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling, 4th. International Heat Transfer Conference, Paris, France, 1970.
↑ Lemmert, M., Chawla, J. M., Influence of flow velocity on surface boiling heat transfer coefficient. Heat Transfer and Boiling, Academic Press, 1977.
↑ Lavieville, J., Quemerais, E., Mimouni, S., Boucker, M., Mechitoua, N., NEPTUNE CFD V1.0 theory manual. NEPTUNE report Nept_2004_L1, 2(3), 2006.
↑ Ranz, W. E., Marshall, W. R. Jr., Evaporation from Drops, Part 2. Chemical Engineering Progress, Vol. 48, No. 4, pp. 173-180, 1952.

Cite error: <ref> tag with name "Incropera2002a" defined in <references> is not used in prior text.
Cite error: <ref> tag with name "Incropera2002b" defined in <references> is not used in prior text.

[RhieChow1983-1] Rhie, C. M., Chow, W. L., Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal 21 (11), pp. 1525-1532, 1983.

[Patankar1980-2] 2.0 ^2.1 Patankar, S. V., Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation. 1980.

[Majumdar1988-3] Majumdar, S., Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids. Numerical Heat Transfer 13 (1), pp. 125-132, 1988.

[Choi1999-4] Choi, S. K., Note on the use of momentum interpolation method for unsteady flows. Numerical Heat Transfer: Part A: Applications 36 (5), pp. 545-550, 1999.

[Pascau2011-5] 5.0 ^5.1 ^5.2 Pascau, A., Cell face velocity alternatives in a structured colocated grid for the unsteady navier-stokes equations. International Journal for Numerical Methods in Fluids 65 (7), pp. 812-833, 2011.

[Kunz1998-6] 6.0 ^6.1 Kunz, R. F., Siebert, B. W., Cope, W. K., Foster, N. F., Antal, S. P., Ettorre, S. M., A coupled phasic exchange algorithm for three–dimensional multi–field analysis of heated flows with mass transfer. Computers & Fluids 27 (7), pp. 741–768, 1998.

[VasquezIvanov2000-7] 7.0 ^7.1 Vasquez, S. A., Ivanov, V. A., A phase coupled method for solving multiphase problems in unstructured meshes. In: Proceedings of ASME FEDSM’00: ASME 2000 fluids engineering division summer meeting, Boston. pp. 743-748, 2000.

[KurulPodowski1991-8] 8.0 ^8.1 Kurul, N., Podowski, M. Z., On the Modeling of Multidimensional Effects in Boiling Channels. In: ANS Proceedings of 1991 National Heat Transfer Conference, Minneapolis, Minnesota, Vol. 5. pp. 30-40, 1991. ISBN: 0-89448-162-2

[DelValleKenning1985-9] 9.0 ^9.1 Del Valle M., V. H., Kenning, D. B. R., Subcooled flow boiling at high heat flux. Int. J. Heat Mass Transfer, Vol. 28, No. 10, pp.1907-1920, 1985.

[McAdams1942-10] McAdams, W. H., Heat Transmission, 2nd Edition, McGraw-Hill, New York, 1942.

[Cole1960-11] Cole, R., A photographic study of pool boiling in the region of CHF. AIChEJ, 6 pp. 533-542, 1960.

[TolubinskiKostanchuk1970-12] Tolubinski, V. I., Kostanchuk, D. M., Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling, 4th. International Heat Transfer Conference, Paris, France, 1970.

[LemmertChawla1977-13] Lemmert, M., Chawla, J. M., Influence of flow velocity on surface boiling heat transfer coefficient. Heat Transfer and Boiling, Academic Press, 1977.

[Lavieville2006-14] Lavieville, J., Quemerais, E., Mimouni, S., Boucker, M., Mechitoua, N., NEPTUNE CFD V1.0 theory manual. NEPTUNE report Nept_2004_L1, 2(3), 2006.

[RanzMarshall1952-15] Ranz, W. E., Marshall, W. R. Jr., Evaporation from Drops, Part 2. Chemical Engineering Progress, Vol. 48, No. 4, pp. 173-180, 1952.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

@@ Line 1,797: / Line 1,797: @@
 <ref name="Kunz1998"> Kunz, R. F., Siebert, B. W., Cope, W. K., Foster, N. F., Antal, S. P., Ettorre, S. M., A coupled phasic exchange algorithm for three–dimensional multi–field analysis of heated flows with mass transfer. Computers & Fluids 27 (7), pp. 741–768, 1998. </ref>
 <ref name="VasquezIvanov2000"> Vasquez, S. A., Ivanov, V. A., A phase coupled method for solving multiphase problems in unstructured meshes. In: Proceedings of ASME FEDSM’00: ASME 2000 fluids engineering division summer meeting, Boston. pp. 743-748, 2000. </ref>
-<ref name="Incropera2002"> Incropera, F. P., DeWitt, D. P., Fundamentals of Heat and Mass Transfer, Chapter 8, p. 470, eq. 8.19. John Wiley & Sons, Inc. 2002. </ref>
+<ref name="Incropera2002a"> Incropera, F. P., DeWitt, D. P., Fundamentals of Heat and Mass Transfer, Chapter 8, p. 470, eq. 8.19. John Wiley & Sons, Inc. 2002. </ref>
-<ref name="Incropera2002"> Incropera, F. P., DeWitt, D. P., Fundamentals of Heat and Mass Transfer, Chapter 8, p. 470, eq. 8.21. John Wiley & Sons, Inc. 2002. </ref>
+<ref name="Incropera2002b"> Incropera, F. P., DeWitt, D. P., Fundamentals of Heat and Mass Transfer, Chapter 8, p. 470, eq. 8.21. John Wiley & Sons, Inc. 2002. </ref>
 <ref name="KurulPodowski1991"> Kurul, N., Podowski, M. Z., On the Modeling of Multidimensional Effects in Boiling Channels. In: ANS Proceedings of 1991 National Heat Transfer Conference, Minneapolis, Minnesota, Vol. 5. pp. 30-40, 1991. ISBN: 0-89448-162-2 </ref>
 <ref name="DelValleKenning1985"> Del Valle M., V. H., Kenning, D. B. R., Subcooled flow boiling at high heat flux. Int. J. Heat Mass Transfer, Vol. 28, No. 10, pp.1907-1920, 1985. </ref>

Difference between revisions of "Kharon theory manual"

Revision as of 11:24, 1 April 2019

Contents

Conservation of Mass

Conservation of Momentum

Conservation of Total Enthalpy

Face Velocity Equation

Phase-Coupled Face Velocity Equation

Pressure & Volume Fraction Correction

Solution Algorithm

Momentum Transfer

Pressure loss

Interface momentum transfer

Heat and Mass Transfer

Surface heat transfer

Interfacial heat transfer

Mass transfer

References

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Tools