Kharon theory manual

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]\, dV = \iiint\limits_V \left(\gamma_{pq} + s_{\text{mass},q}\right)\, dV$

(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]\, dA = \iiint\limits_V \left(\gamma_{pq} + s_{\text{mass},q}\right)\, dV$

(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/m3s],
$\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/m3],
$A$	face area [m2],
$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/m3s],
$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]\, dV = \iiint\limits_V \left[\text{max}\left(\gamma_{pq},0\right)u_p - \text{max}\left(-\gamma_{pq},0\right)u_q\right]\, dV$

(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)\, dV$

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]\, dA = \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]\, dA$

(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]\, dA$

(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]\, dA$

(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]\, dA$

(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.

Finally, denoting the central coefficient by $a_P = 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)$ 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} + K_{\text{c1}} & -K_{\text{c1}} & -a_{\text{e,c1},1} & 0 & \cdots & 0 \\ -K_{\text{c1}} & a_{\text{P,c1},2} + K_{\text{c1}} & 0 & -a_{\text{e,c1},2} & \ddots & \vdots \\ -a_{\text{w,c2},1} & 0 & a_{\text{P,c2},1} + K_{\text{c2}} & -K_{\text{c2}} & -a_{\text{e,c2},1} & 0 \\ 0 & -a_{\text{w,c2},2} & -K_{\text{c2}} & a_{\text{P,c2},2} + K_{\text{c2}} & 0 & -a_{\text{e,c2},2} \\ \vdots & \ddots & -a_{\text{w,c3},1} & 0 & a_{\text{P,c3},1} + K_{\text{c3}} & -K_{\text{c3}} \\ 0 & \cdots & 0 & -a_{\text{w,c3},2} & -K_{\text{c3}} & a_{\text{P,c3},2} + K_{\text{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}$

(20)

where:

$S_{\text{mom},\text{c1},q} = \Big\{\big[\text{max}\left(-\Gamma_{pq},0\right) + \text{max}\left(-S_{\text{mass},q},0\right)\big]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}\Big\}_{\text{c1}}$

and

$\alpha_q$		volume fraction of phase q [-],
$\gamma_{pq}$		volumetric mass transfer rate from phase p to q [kg/m3s],
$\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/m2],
$\rho_q$		density of phase q [kg/m3],
$A$		face area [m2],
$C$		inertial loss coefficient [1/m],
$F$		face mass flow rate [kg/s],
$g$		acceleration of gravity [m/s2],
$k$		volumetric interphase drag coefficient [kg/m3s],
$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/m3s],
$S_{\text{mass},q}$		mass source to phase q [kg/s],
$s_{\text{mom},q}$		volumetric momentum source to phase q [N/m3],
$S_{\text{mom},q}$		momentum source to phase q [N],
$u_q$		velocity of phase q [m/s],
$V$		cell volume [m3],
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)$

(21)

$= \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)$

(22)

$= \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}$

(23)

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}$

(24)

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 23 and Equation 24 to Equation 22 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)$

(25)

$= \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)$

(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}$

Equation 26 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]\, dV$

(27)

$= \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\}\, dV$

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]\, dA$

(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\}\, dV$

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]\, dA$

(29)

$= \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}$

(30)

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]\, dA$

(31)

$= \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 31 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]\, dA$

(32)

$= \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}$

(33)

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]\, dA$

(34)

$= \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 34 into Equation 28 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}$

(35)

$= \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\}\, dV$

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] dV$

(36)

$= \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] dV$

(37)

$= -\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] dV = - \left( F_{\text{D}} + F_{\text{IPD}} \right) u_{\text{P},q}$

(38)

where the drag force, $F_{\text{D}}$, and the interphase drag force, $F_{\text{IPD}}$, 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$

(39)

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

(40)

Volume integral of volumetric enthalpy source to phase q is simply

$\iiint\limits_V s_{h,q} dV = S_{h,q}$

(41)

Combining equations Equation 35, Equation 36, Equation 37, Equation 38 and Equation 41 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}$

(42)

$= \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{IPD}} \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}$

(43)

$= -\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{IPD}} \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 43

$\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}$

(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{IPD}} \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}$

(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{IPD}} \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}$

(46)

$= \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{IPD}} \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/m3s],
$\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/m3],
$A$		face area [m2],
$C$		inertial loss coefficient [1/m],
$F$		face mass flow rate [kg/s],
$F_{\text{D}}$		drag force [N],
$F_{\text{IPD}}$		interphase drag force [N],
$\mathbf{g}$		acceleration of gravity (vector) [m/s2],
$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/m3],
$S_{h,q}$		energy source to phase q [W],
$s_{\text{mass},q}$		volumetric mass source to phase q [kg/m3s],
$S_{\text{mass},q}$		mass source to phase q [kg/s],
$u_q$		velocity of phase q [m/s],
$V$		cell volume [m3],
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 46 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 remedies, at least the most widely known, was proposed by Rhie & Chow. However, they did not provide a proper derivation for their formulation. The following aims to provide such a derivation. In addition, although not directly evident, the results are 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}}$

(47)

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}^*$

(48)

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}^*$

(49)

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}^*$

(50)

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

$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}^*$

(51)

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}^*}$

(52)

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}$

(53)

Solving $\frac{\overline{\sum_{nb=1}^N a_{\text{nb}} u_{\text{nb},q} + S_{\text{mom},f,q}}} {a_f}$ from Equation 52 and substituting it to Equation 51 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( \overline{u_{f,q}^*} - u_{f,q}^* \right)$

(54)