Difference between revisions of "Kharon theory manual"

Revision as of 13:33, 2 January 2019

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

(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/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{IPD}}$		interphase 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 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, 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}}$

(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( u_{f,q}^* - \overline{u_{f,q}^*} \right)$

(54)

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

(55)

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

(56)

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

(57)

Casting Equation 57 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}$

(58)

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

(59)

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

(60)

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

(61)

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

(62)

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

(63)

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

(64)

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

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

(65)

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

(66)

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

(67)

Solving $\overline{\left( \mathbf{A} - \mathbf{R} \right)^{-1} \overrightarrow{H_f}}$ from Equation 66 and substituting it to Equation 65 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)$

(68)

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

(69)

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

(70)

Linear interpolation of the first term inside the brackets, in Equation 69, 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}$

(71)

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

(72)

Substituting the results of Equation 71 and Equation 72 to Equation 69 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)$

(73)

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

(74)

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

(75)

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

(76)

The relation between pressure and velocity corrections is obtained from equation Equation 73 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}$

(77)

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

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

(78)

Introducing the corrected face mass flow rates from Equation 74 into Equation 78.

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

(79)

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

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

(80)

Substituting Equation 77 into Equation 80 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}^*

(81)

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

(82)

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

(83)

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 82

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

(84)

Writing Equation 84 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}$

(85)

A single system of equations is formed from Equation 85 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}$

(86)

For some linear solvers it might be beneficial to load the equations in a different order. The ordering shown in Equation 20 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 86 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)$

(87)

phase volume fractions are corrected with

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

(88)

and pressure is corrected with

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

(89)

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

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

Wall boiling model

The RPI wall boiling model is a mechanistic model developed at Rensselaer Polytechnic Institute [8]. 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 the wall heat flux, $q_{\text{w}}$ [W/m²], is partitioned into three parts, as follows

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

(90)

where:

$q_{\text{c}}$		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²].

Next, an additional partitioning is introduced to extend applicability towards high void fractions:

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

(91)

Interfacial heat transfer model

References

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

[2] Patankar, S. V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation.

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

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

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

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

[7] Vasquez, S. A., Ivanov, V. A., 2000. 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.

[8] Kurul, N., Podowski, M. Z., 1991. 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. ISBN: 0-89448-162-2

@@ Line 1,318: / Line 1,318: @@
 {{NumBlk|:|
 <math>
-q_{\text{w}} = \left( q_{\text{c}} + q_{\text{q}} + q_{\text{e}} \right) + \left( 1 - \right)
+q_{\text{w}} = \left( q_{\text{c}} + q_{\text{q}} + q_{\text{e}} \right) + \left( 1 - f_{\text{liquid}} \right)
 </math>
 |<xr id="eqn:WallHeatFluxPartitioning2" nolink />}}

Difference between revisions of "Kharon theory manual"

Revision as of 13:33, 2 January 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

Heat Transfer

Wall boiling model

Interfacial heat transfer model

References

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Tools