---
title: "Pressure-Temperature-Equilibrium Multi-Material Algebraic-VOF Framework"
---

# Pressure-Temperature-Equilibrium Multi-Material Algebraic-VOF Framework

This chapter specifies a pressure-temperature-equilibrium (PTE),
single-velocity model for compressible multi-material, multispecies flow. Every
material in a cell shares one pressure, one temperature, and one velocity, but
retains its own equation of state (EOS) and material composition. Materials
are immiscible and are tracked by the algebraic volume-of-fluid (VOF) method.
Passive species diffusion occurs only among species belonging to the same
material.

The conserved species densities, mixture momentum, and mixture total
enthalpy are advanced by their transport equations. Each material mass density
is the sum of its species densities, and mixture density is the sum over all
materials. Pressure correction enforces the material-volume constraint
$\sum_m\alpha_m=1$.

The model uses the following assumptions:

1. One common pressure, temperature, and velocity in each cell. Material
   closures are evaluated at this common $(p,T)$ and with their own
   material compositions.
2. Material interfaces are tracked by VOF fields. Passive species diffusion
   does not transfer species between different materials.
3. Species belong to exactly one material set and mix only with species owned
   by that same material.
4. Each material has one EOS closure for its internally miscible species.
   Mixture properties are assembled from the independently evaluated material
   states and their material masses.
5. Chemistry is a finite-rate conservative source closure. Phase change is a
   conservative inter-material transfer. A resolved surface rate is determined
   by the interfacial heat balance, while a volumetric cavitation rate is
   supplied by a separate constitutive closure. Common pressure and temperature
   do not imply chemical-potential or phase equilibrium.

Every species balance is solved, including the sole species of a one-species
material. Summing the species balances within a material gives its material
mass balance; summing over all materials gives mixture continuity. These
summed balances are consequences of the species equations, not additional
equations to solve. Material volume fractions are derived from the species
masses and the material EOS densities.

The discretization uses a fixed control-volume mesh, the first-order backward
differentiation formula (BDF1), and a segregated pressure-based iteration. The
thermal equation advances total enthalpy; the static-enthalpy equation is its
equivalent continuum form. Material occupancy may use limited second-order
upwind (SOU) reconstruction or the compressive MIG reconstruction without a
separate artificial compression velocity. SOU changes the spatial face
approximation; BDF1 remains first-order in time.

Molecular transport uses volume-fraction-weighted viscosity and thermal
conductivity, Fourier heat conduction, and mixture-averaged species diffusion
in each material. Capillarity uses constant pairwise surface tension, CELESTE
interface curvature, and balanced pressure--surface-force coupling. Surface
phase change uses the Kunkelmann interfacial heat-balance method and the
conservative transfer relations defined below. Pressure-driven volumetric
transfer is treated as a distinct local rate closure.


## Indices, Sets, and Ownership Maps

This framework distinguishes quantities defined at the material level from
quantities defined at the global mixture level. A material is a collection of
species tracked by one VOF field. The indices used below are:

1. $m = 1,...,N_m$ (index of materials)
2. $k = 1,...,N_s$ (index of species), with
   $N_{s,m}=|\mathcal K_m|$ species owned by material $m$
3. $\Omega$ is the control-volume size of the cell

Species ownership sets are disjoint. The sets $\mathcal{K}_m$ contain the
species indices owned by material $m$. Every species belongs to exactly one
material, with no overlap and no unowned species:

```{math}
:label: chap13_material_species_sets

\mathcal{K}_m\subset\{1,\dots,N_s\},
\qquad
\mathcal{K}_m\cap\mathcal{K}_n=\emptyset\ (m\neq n),
\qquad
\bigcup_{m=1}^{N_m}\mathcal{K}_m=\{1,\dots,N_s\}
```

The volume fraction $\alpha_m$ is the fraction of the control volume occupied
by material $m$. Each fraction lies between zero and one, and their sum is one:

```{math}
:label: chap13_alpha_constraints

0\le\alpha_m\le1,
\qquad
\sum_{m=1}^{N_m}\alpha_m=1
```

The owner-material map $m(k)$ is defined by
$k\in\mathcal{K}_{m(k)}$. It identifies the material that owns species $k$
and hence the composition vector containing that species.

All $N_{s,m}$ species balances are retained and solved, including the sole
species balance of a one-species material.

The cell-integrated material mass ($M_m$), species mass ($M_k$), and their
conservative densities are

```{math}
:label: chap13_cell_masses

\begin{aligned}
z_{mk}&=\frac{M_k}{\Omega}=\alpha_m\rho_m\widetilde Y_k^{(m)},
\qquad k\in\mathcal K_m,\\
q_m&=\frac{M_m}{\Omega}=\alpha_m\rho_m.
\end{aligned}
```

Here $\rho_m$ is the intrinsic density of material $m$, and
$\widetilde{Y}^{(m)}_k$ is the mass fraction of species $k$ within its owner
material. The transported quantity is $z_{mk}$. Its material sum defines the
partial mass density at every time and iteration level:

```{math}
:label: chap13_material_species_mass_sum

q_m=\sum_{k\in\mathcal K_m}z_{mk}.
```

For a one-species material, $z_{mk}=q_m$ and
$\widetilde Y_k^{(m)}=1$ wherever the material is present. In general, the
species mass fraction within material $m$ is

```{math}
:label: chap13_local_mass_fraction_definition

\widetilde{Y}^{(m)}_k=\frac{M_k}{M_m},
\qquad
k\in\mathcal{K}_m,
\qquad
M_m>0
```

The vector $\widetilde{\vec{Y}}^{(m)}$ is called the **composition of
material $m$**. It contains the species mass fractions within that material
and has its own unit-sum constraint. The species fractions sum to one
separately within each material. Material composition does
not describe how much of each material occupies the cell; that distribution
is represented by the volume fractions $\alpha_m$. The related material mass
fractions are $w_m$. The composition is conditional on material $m$ being
present and is physically defined only where $q_m>0$. The mixture species
fractions $Y_k$ defined below are derived from the material masses and material
compositions; they are not the composition supplied to an individual material
EOS. For $k\notin\mathcal{K}_m$, the quantity is not part of material $m$'s
thermodynamic state. Define the active-material set

```{math}
:label: chap13_active_material_set

\mathcal A(\mathbf x,t)=\{m:q_m(\mathbf x,t)>0\}.
```

Every physical mixture sum and conditional transport operation involving a
material composition or evaluated material property is taken over
$\mathcal A$, even when an all-material sum is written for compactness. An
absent term is omitted before its conditional property is evaluated. A
generic exact-zero cell has no prospective EOS state; a material introduced by
transport joins $\mathcal A$ only after its species predictor has supplied
a complete positive state. A physical creation model must likewise supply a
complete scoped source state rather than borrowing dormant cell properties.
The derived mixture species mass fraction is

```{math}
:label: chap13_global_from_internal

Y_k=\frac{M_k}{\sum_{n=1}^{N_m}M_n} = \frac{\alpha_{m(k)}\rho_{m(k)}}{\rho}\widetilde{Y}^{(m(k))}_k
```

The mass fraction of material $m$ in the full mixture is

```{math}
:label: chap13_material_mass_fraction

w_m=\frac{\alpha_m\rho_m}{\rho},
\qquad
\sum_{m=1}^{N_m}w_m=1
```

The species mass fraction within its material maps to the mixture species
mass fraction as

```{math}
:label: chap13_local_to_global_transformation

Y_k=w_{m(k)}\widetilde{Y}^{(m(k))}_k
```

The species fractions belonging to one material sum to its mixture mass
fraction:

```{math}
:label: chap13_vof_species_manifold

\sum_{k\in\mathcal K_m}Y_k
=w_m\sum_{k\in\mathcal K_m}\widetilde Y_k^{(m)}
=w_m.
```

This mapping is used for mixture reporting and for identities obtained by
summing the material states. It introduces no additional transported or
history variable. The mixture species vector is not used to evaluate species
diffusion in a material, reconstructed to obtain a material face composition,
or supplied to a material EOS.

For an exact admissible state, these definitions have the properties:

```{math}
:label: chap13_species_fraction_properties
Y_k\ge0,
\qquad\sum_{k=1}^{N_s}Y_k=1,
\qquad\widetilde{Y}^{(m)}_k\ge0,
\qquad\sum_{k\in\mathcal{K}_m}\widetilde{Y}^{(m)}_k=1
```

:::{admonition} Meaning of simplex in this chapter

For a collection of $N$ fractions, the **simplex** is the admissible set

```{math}
\mathcal S_N
=
\left\{
\mathbf x\in\mathbb R^N:
x_i\ge0,
\ \sum_{i=1}^{N}x_i=1
\right\}.
```

In plain language, every fraction is nonnegative and all the fractions add to
one. The **material-volume simplex** applies this requirement to
$(\alpha_1,\ldots,\alpha_{N_m})$. Each active material has a separate
**composition simplex** for the species fractions
$\widetilde{\vec{Y}}^{(m)}$ that belong to that material. Species fractions
from different materials are not combined into one composition simplex.

When only the sum-to-one property is used in a derivation, that property is
stated explicitly. Saying that a numerical operation *preserves the simplex*
means that it preserves both nonnegativity and the unit sum.
:::

The amount of material is carried by $q_m$, not by its conditional
composition. At $q_m=0$, any extension of
$\widetilde{\vec{Y}}^{(m)}$ is arbitrary. A conditional-state SOU increment
is therefore used only when material $m$ is present throughout the donor's
gradient stencil. Otherwise that material uses the upwind-cell composition.
This support decision is made independently for every material: an absent
material does not reduce the reconstruction order of another material with a
fully supported stencil. Advective composition is used only when the material
face mass flux is nonzero, and molecular diffusion is evaluated only across a
face supported by that material on both sides.


## State Variables and Closure Fields

Each material owns its species set, equation of state, and transport
properties. Three roles must remain distinct: conserved histories, variables
used by the segregated solves, and thermodynamic closure fields.


### Conserved quantities, segregated unknowns, and closure

The material-species equations use $z_{ms}$ as the segregated linear-system
unknowns and as accepted time histories. Their material sums define $q_m$ at
each time and iteration level. The material balance follows by summing its
species balances. Material volume fraction and conditional composition are
derived closure views. Momentum and total enthalpy retain the primitive-like
segregated unknowns shown below.

For each present material, the species densities are nonnegative and their
sum is positive. The composition satisfies

```{math}
:label: chap13_material_internal_species

\widetilde{Y}^{(m)}_s\ge0,
\qquad
\sum_{s\in\mathcal{K}_m}\widetilde{Y}^{(m)}_s=1.
```

It is obtained directly from the species densities:

```{math}
:label: chap13_material_query_composition

\widetilde Y_s^{(m)}=\frac{z_{ms}}{q_m},
\qquad
q_m=\sum_{s\in\mathcal K_m}z_{ms}>0.
```

This relation expresses the mass fractions within a material; it does not
change any species mass. If $q_m=0$, all its species densities are zero and
the conditional composition is undefined.

| Equation | Conserved state/history | Segregated unknown | Lagged coefficients |
| --- | --- | --- | --- |
| Material species | $z_{ms}$ | every $z_{ms}^{*,k}$ | volume flux, material and composition reconstruction, diffusion |
| Momentum | $\rho\vec{u}$ | $\vec{u}^{*,k}$ | density, face mass flux, transport properties |
| Total energy | $\rho E=\rho h_t-p$ | $h_t^{*,k}$ | thermodynamic state, species and mixture face fluxes |

The quantities $q_m=\sum_s z_{ms}$, $\rho=\sum_m q_m$,
$\widetilde Y_s^{(m)}=z_{ms}/q_m$, and $\alpha_m=q_m/\rho_m$ are derived
at each iteration. They are not additional segregated unknowns.

Common $p,T$ and each material EOS close these balances. The resulting
closure fields include $\rho_m,h_m,e_m$ and their derivatives, followed by
$\rho,w_m$, and the derived mixture species fractions $Y_k$.

(chap13_outer_iteration)=
#### Outer iteration

The general pressure-based iteration is defined in
{numref}`Chapter %s <chap_solution_algorithm>`. For CMIM, the accepted state at
time level $n$ supplies the time histories and the initial iterate at the new
time. A superscript $k$ denotes the incoming outer iterate, $*$ denotes the
result of the owning predictor, and a prime denotes a pressure correction. A
$\mathrm{pred},k$ label denotes the thermodynamic state evaluated after the
predictors at the incoming pressure $p^k$; it is not another predictor solve.
When both indices occur, $s\in\mathcal K_m$ denotes a species and the
superscript $k$ denotes the outer iteration.

Each segregated predictor has the finite-volume row form established in
{numref}`Chapter %s <chap_discretized_equations_unstructured>`:

```{math}
{}^{(\phi)}\!a_P^k\phi_P^{*,k}
=
\sum_{N\in\mathcal N(P)}{}^{(\phi)}\!a_{P,N}^k\phi_N^{*,k}
+{}^{(\phi)}\!b_P^k,
\qquad
\phi\in\{z_{ms},u_i,h_t\}.
```

The coefficients and right-hand side are assembled from the accepted time
histories and the same incoming iterate $k$, then held fixed during the linear
solve. One predictor does not use another predictor's newly solved value. The
outer iteration supplies their coupling. The experimental sequential energy
update in {ref}`chap13_material_energy_coupling` changes this dependency; it
is separate from the standard sequence described here.

The information flow through one CMIM outer iteration is

```text
Accepted species, momentum, and energy histories at time n
                         |
                         v
Incoming outer iterate k and its derived material states
                         |
                         v
Independent predictors, with coefficients from iterate k
  species: z_ms -> z_ms*; derive q_m* = sum_s z_ms*
  momentum: u -> u*, predicted face volume flow
  total enthalpy: h_t -> h_t*
                         |
                         v
Pre-pressure closure at p^k with the predictor results
  solve T; evaluate rho_m; derive alpha_m = q_m*/rho_m
                         |
                         v
Pressure correction, followed by corrections to u and face volume flow
  hold z_ms* and h_t* fixed
                         |
                         v
Post-pressure closure
  solve T; evaluate rho_m; derive alpha_m = q_m*/rho_m
                         |
                         v
Nonlinear convergence and admissibility checks
  satisfied: accept time n+1
  otherwise: repeat with the new outer iterate
```

Only a converged admissible state becomes time level $n+1$. The complete
sequence is given in {ref}`chap13_segregated_nonlinear_algorithm`.

(chap13_material_absence_and_appearance)=
### Material absence and appearance

Material composition and the material EOS state are conditional on material
mass. In an admissible cell, material $m$ is absent exactly when its complete
conservative block is zero:

```{math}
q_{m,P}=0,
\qquad
z_{ms,P}=0
\quad(s\in\mathcal K_m).
```

At this state, the material composition and EOS state are undefined.
Active-material sums include materials with $q_{m,P}>0$, and EOS evaluations
are applied to those materials. Every finite positive $q_{m,P}$ belongs to the
active conservative state and has the complete species block specified above.

Floating-point transport can leave negligible positive material densities in
otherwise single-material cells. Before publishing the species predictor, CMIM
zeros an entire material block when its total density is below a configurable
cutoff, defaulting to $10^{-20}$ kg/m$^3$. Individual trace species within a
retained material are preserved. The same policy applies after validating
current and previous restart states. Consequently, the exact-zero test above
identifies intentional material absence before any EOS query.

This cleanup deliberately removes mass. For effective cutoff $q_c$, the mass removed
from an affected cell of volume $V_P$ is less than $q_c V_P$ per material and
cleanup operation. Velocity and specific total enthalpy are retained; their
extensive conserved quantities therefore change with the removed mass.
The cutoff must be small relative to the case's physical mass scales, and
verification must assess its conservation effect. Setting the input cutoff to
zero retains only the existing floating-point underflow protection.

The species balances determine material appearance. Their cell rows combine accepted history, material-bearing face fluxes, and physical
sources. Starting from an accepted exact-zero block, an inflow or a
material-creation source can produce positive predicted material mass. Matched
material and species face fluxes carry the donor composition into the receiving
cell, and the cell balance combines contributions from multiple faces. A
physical source that creates material supplies its material mass, every owned
species contribution, the corresponding energy contribution, and any momentum
required by the source model. The material/species source identity is given in
Eq. {eq}`chap13_material_transfer_source_constraint`.

After the predictors, a positive material sum
$q_{m,P}^{*,k}=\sum_{s\in\mathcal K_m}z_{ms,P}^{*,k}>0$ defines the material
composition through Eq. {eq}`chap13_material_query_composition`.
At the pre-pressure closure, evaluation at
$(p_P^k,T_P^{\mathrm{pred},k})$ supplies $\rho_{m,P}^{\mathrm{pred},k}$, and
$\alpha_{m,P}^{\mathrm{pred},k}=q_{m,P}^{*,k}/
\rho_{m,P}^{\mathrm{pred},k}$ is derived for the pressure equation. After
pressure correction, the common-temperature closure rebuilds
$\rho_{m,P}^{k+1}$ and $\alpha_{m,P}^{k+1}$. Neither density nor volume
fraction is an independently solved predictor. Detailed face configurations,
reconstruction support, and pressure-row participation are given in
{ref}`chap13_material_appearance_face_reference`.

### Derived Fields

| Quantity | Symbol | Definition |
| --- | --- | --- |
| Material mass density | $q_m$ | $q_m=\sum_{s\in\mathcal K_m}z_{ms}$ |
| Material volume fraction | $\alpha_m$ | $\alpha_m=q_m/\rho_m$ for a present material; zero otherwise |
| Partial specific enthalpy | $\bar h_{mk}$ | $\bar h_{mk}=\bar h_{mk}(\mathbf{s}_m)$ for $k\in\mathcal{K}_m$ |
| Material density | $\rho_m$ | Eq. {eq}`chap13_material_density` |
| Material enthalpy | $h_m$ | Eq. {eq}`chap13_material_enthalpy` |
| Material internal energy | $e_m$ | Eq. {eq}`chap13_material_internal_energy` |
| Mixture density | $\rho$ | Eq. {eq}`chap13_mixture_density` |
| Mixture static thermodynamic enthalpy | $h$ | Eq. {eq}`chap13_mixture_enthalpy` |
| Mixture internal energy | $e$ | Eq. {eq}`chap13_mixture_internal_energy` |
| Mixture total energy | $E$ | Eq. {eq}`chap13_mixture_total_energy` |
| Mixture species mass fraction | $Y_k$ | Eq. {eq}`chap13_global_from_internal` |
| Material mass fraction | $w_m$ | Eq. {eq}`chap13_material_mass_fraction` |
| Fixed-enthalpy volume compressibility | $\chi$ | Eq. {eq}`chap13_fixed_q_volume_response` |


## Thermodynamic Definitions

Material thermodynamic properties are evaluated from pressure, temperature,
and material composition. The materials share the same cell pressure and
temperature, but each material EOS receives only its own composition.

It is useful to distinguish a material constitutive model from an evaluated
material state. Let $\mathcal{E}_m$ denote the EOS mapping for material $m$.
A cell-local thermodynamic state is the result of evaluating that mapping:

```{math}
:label: chap13_eos_evaluator_state

\mathbf{s}_m(\mathbf{x},t)
=
\mathcal{E}_m\!\left(
p(\mathbf{x},t),
T(\mathbf{x},t),
\widetilde{\vec{Y}}^{(m)}(\mathbf{x},t)
\right).
```

For every active material, the state $\mathbf{s}_m$ supplies finite
$\rho_m>0$, $h_m$, $e_m$, the partial specific enthalpies required by species
diffusion and phase transfer, the physical sound speed $c_m$ at fixed material
composition, and the thermodynamic derivatives required by the PTE closure.
It also reports whether
the state lies on the material's assigned phase branch and inside its valid
thermodynamic domain. The constitutive state is determined only by the local
$(p,T,\widetilde{\vec{Y}}^{(m)})$ and the material model $\mathcal E_m$.

The evaluated state supplies the intrinsic material density and static
thermodynamic enthalpy:

```{math}
:label: chap13_material_density

\rho_m
=
\rho_m\!\left(
p,T,\widetilde{\vec{Y}}^{(m)};\mathcal{E}_m
\right),
```

```{math}
:label: chap13_material_enthalpy

h_m
=
h_m\!\left(
p,T,\widetilde{\vec{Y}}^{(m)};\mathcal{E}_m
\right).
```

It also supplies the fixed-composition derivatives

```{math}
:label: chap13_cmim_material_eos_partial_derivatives

\rho_{m,p}
=
\left(\frac{\partial\rho_m}{\partial p}\right)_{T,\widetilde{\vec{Y}}^{(m)}},
\qquad
\rho_{m,T}
=
\left(\frac{\partial\rho_m}{\partial T}\right)_{p,\widetilde{\vec{Y}}^{(m)}},
\qquad
h_{m,p}
=
\left(\frac{\partial h_m}{\partial p}\right)_{T,\widetilde{\vec{Y}}^{(m)}},
\qquad
h_{m,T}
=
\left(\frac{\partial h_m}{\partial T}\right)_{p,\widetilde{\vec{Y}}^{(m)}}.
```

The mapping $\mathcal E_m$ returns these derivatives with $\rho_m$ and $h_m$,
evaluated at the same $(p,T,\widetilde{\vec{Y}}^{(m)})$ on the assigned EOS
branch.

For the species-diffusion energy flux, let
$H_m=M_mh_m$ be the extensive material enthalpy. The same evaluated state
supplies the partial specific enthalpy

```{math}
:label: chap13_species_eos

\bar h_{mk}
=
\left(
\frac{\partial H_m}{\partial M_k}
\right)_{p,T,M_{\ell\ne k},\,\ell\in\mathcal K_m}
\qquad k\in\mathcal{K}_m.
```

Extensivity gives

```{math}
:label: chap13_partial_enthalpy_extensivity

H_m
=
\sum_{k\in\mathcal K_m}M_k\bar h_{mk},
\qquad
h_m
=
\sum_{k\in\mathcal K_m}\widetilde Y_k^{(m)}\bar h_{mk}.
```

The partial specific enthalpy is evaluated at the complete material
composition with the same reference convention as $h_m$.

Each $\mathcal E_m$ defines the mixing rule for the species inside material $m$.
There is no cross-material EOS and no second reconstruction of a material
state from mixture species fractions. The mixture volume is the sum of the
independently evaluated material volumes:

```{math}
:label: chap13_geometric_material_volume

\Omega=\sum_m\frac{M_m}{\rho_m},
\qquad
\alpha_m=\frac{q_m}{\rho_m},
\qquad
\sum_m\frac{q_m}{\rho_m}=1.
```

The material internal energy is related to its static enthalpy by

```{math}
:label: chap13_material_internal_energy

e_m=h_m-\frac{p}{\rho_m}.
```

The mixture level is the level of the common pressure, velocity, and
temperature fields. Each material nevertheless retains its own evaluated EOS
state. The mixture density is mass per mixture volume, so each material
contributes according to the volume fraction it occupies:

```{math}
:label: chap13_mixture_density

\rho=\sum_{m=1}^{N_m}q_m=\sum_{m=1}^{N_m}\alpha_m\rho_m
```

The mixture static thermodynamic enthalpy is a specific quantity per unit
mixture mass, so it is averaged with the material mass fractions:

```{math}
:label: chap13_mixture_enthalpy

h=\sum_{m=1}^{N_m}\frac{\alpha_m \rho_m}{\rho} h_m = \sum_{m=1}^{N_m} w_m h_m
```

The mixture internal-energy density follows from the same material masses:

```{math}
:label: chap13_mixture_internal_energy

\rho e
=
\sum_{m=1}^{N_m}\alpha_m\rho_m e_m
=
\rho h-p,
```

where the final equality uses $\sum_m\alpha_m=1$. Defining
$K=\frac12\lVert\vec{u}\rVert^2$, the mixture total energy is

```{math}
:label: chap13_mixture_total_energy

E=e+K,
\qquad
\rho h_t-p=\rho E.
```

The total enthalpy is defined as the static thermodynamic enthalpy plus kinetic
energy:

```{math}
:label: chap13_total_enthalpy_def

h_t=h+\frac{1}{2}\lVert\vec{u}\rVert^2
```

### Pressure-temperature compatibility

At a fixed transported state, the species densities $z_{ms}$, their material
sums $q_m$, velocity $\vec{u}$, and total enthalpy $h_t$ are known.
For each material in the active set $\mathcal A$, the
$z_{ms}$ block determines $\widetilde{\vec{Y}}^{(m)}$ through
Eq. {eq}`chap13_material_query_composition`. The material EOS branches are
also held fixed. The common pressure and temperature must satisfy the volume
and enthalpy compatibility conditions below.

With the kinetic energy defined above, the mixture density and target static
enthalpy are

```{math}
\rho=\sum_{m\in\mathcal A}q_m,
\qquad
h^{\mathrm{tar}}=h_t-K.
```

A trial $(p,T)$ supplies $\rho_m(p,T,\widetilde{\vec{Y}}^{(m)})$ and
$h_m(p,T,\widetilde{\vec{Y}}^{(m)})$ from each active material EOS. Because
$q_m\Omega$ is the material mass in a cell, $q_m/\rho_m$ is the fraction of
the cell volume occupied by that material. Additive volume requires these
fractions to sum to one. Equation {eq}`chap13_mixture_enthalpy` also requires
$\sum_{m\in\mathcal A}q_mh_m=\rho h^{\mathrm{tar}}$. These two requirements
define the closure residuals

```{math}
:label: chap13_pte_closure_residuals

F_V(p,T)
=
\sum_{m\in\mathcal A}
\frac{q_m}{\rho_m(p,T,\widetilde{\vec{Y}}^{(m)})}-1
=
\sum_{m\in\mathcal A}\alpha_m(p,T)-1
=0,
\qquad
F_H(p,T)
=
\sum_{m\in\mathcal A}q_m
\left[
h_m(p,T,\widetilde{\vec{Y}}^{(m)})-h^{\mathrm{tar}}
\right]
=0.
```

Here $\alpha_m(p,T)=q_m/\rho_m(p,T,\widetilde{\vec{Y}}^{(m)})$ is a trial
EOS-derived volume fraction, not an independent input. The residual $F_V$ is
dimensionless. The residual $F_H$ has units $\mathrm{J/m^3}$, equivalently
$\mathrm{Pa}$; dividing it by $\rho$ gives the corresponding
specific-enthalpy residual in $\mathrm{J/kg}$. These are compatibility
conditions on the transported state, not additional conservation equations.

The same residuals enforce the total-energy representation. Using
$e_m=h_m-p/\rho_m$ gives

```{math}
\left(\sum_{m\in\mathcal A}q_m e_m+\rho K\right)
-\left(\rho h_t-p\right)
=
F_H-pF_V.
```

Therefore, $F_V=F_H=0$ makes the total-energy density reconstructed from the
material EOS states equal to the transported total-energy density while the
material volumes fill the cell. At the coupled root, this enthalpy form is
equivalent to the specific-volume and internal-energy PTE closure discussed in
References [5] and [6].

Assume every active EOS branch is continuous over the common admissible
interval and satisfies $c_{p,m}>0$ and $\rho_{m,p}>0$ there. Then

```{math}
\left(\frac{\partial F_H}{\partial T}\right)_p
=
\sum_{m\in\mathcal A}q_m c_{p,m}>0,
\qquad
\left(\frac{\partial F_V}{\partial p}\right)_T
=
-\sum_{m\in\mathcal A}
\frac{q_m\rho_{m,p}}{\rho_m^2}
<0.
```

Thus, at fixed pressure, the enthalpy residual has at most one temperature
root. A root exists in the selected interval when the active EOS temperature
domains overlap and the endpoint residuals bracket zero. Both residuals
generally depend on both $p$ and $T$.

For a fixed pressure $p$, the common-temperature inversion is the scalar problem

```{math}
:label: chap13_temperature_from_enthalpy

f_H(\theta)
\equiv
\frac{1}{\rho}
\sum_{m\in\mathcal A}q_m
\left[
h_m(p,\theta,\widetilde{\vec{Y}}^{(m)})-h^{\mathrm{tar}}
\right],
\qquad
f_H(T)=0.
```

During this local solve, $p$, $q_m$, $z_{ms}$, $h^{\mathrm{tar}}$, the active
set, and the material EOS branches remain fixed. Each active material EOS is
evaluated at every trial temperature $\theta$.

CMIM uses this inversion both before and after pressure correction. The
pre-pressure evaluation supplies the material densities and derivatives in
Eq. {eq}`chap13_fixed_q_volume_response`. The pressure equation uses the local
derivative of $F_V$ along $F_H=0$ together with the pressure-induced
face-volume-flux response; it is not a local two-variable closure solve. After
pressure correction, the corrected velocity defines $h^{\mathrm{tar}}$, the
scalar inversion supplies the common temperature, and the material EOS states
are evaluated at the corrected $(p,T)$. The complete stage order is given in
{ref}`chap13_segregated_nonlinear_algorithm`.

## Mixture Viscosity and Heat Conduction

For each material with positive partial density $q_m$, its transport model supplies the dynamic viscosity $\mu_m$ and thermal
conductivity $\lambda_m$ at the common cell pressure and temperature and the material-local composition. The cell mixture
properties are the normalized volume-fraction averages

```{math}
\alpha_\Sigma=\sum_m\alpha_m,
\qquad
\mu_{\mathrm{mix}}=\frac{1}{\alpha_\Sigma}\sum_m\alpha_m\mu_m,
\qquad
\lambda_{\mathrm{mix}}=\frac{1}{\alpha_\Sigma}\sum_m\alpha_m\lambda_m.
```

An absent material has $\alpha_m=0$ and does not contribute. At closure, $\alpha_\Sigma=1$; the normalized form is retained during
the nonlinear iteration when the volume-fraction sum may differ from unity.

For turbulent flow, the momentum and thermal operators use the effective mixture coefficients

```{math}
\mu_{\mathrm{eff}}=\mu_{\mathrm{mix}}+\mu_t,
\qquad
\lambda_{\mathrm{eff}}
=\lambda_{\mathrm{mix}}+\frac{\mu_t c_p}{Pr_t},
```

where $\mu_t$ is the eddy viscosity and $Pr_t$ is the turbulent Prandtl number. For laminar flow, $\mu_t=0$. Below, $\mu$ and
$\lambda$ denote these effective coefficients. The turbulent thermal contribution does not add turbulent species diffusion; the
CMIM species-diffusion model remains molecular.

Internal-face momentum viscosity is the arithmetic average of the effective cell viscosities. Fourier conduction instead uses a
distance-weighted harmonic coefficient. For a face $f$ between cells $L$ and $R$, define

```{math}
\theta_f=\frac{d_{Lf}}{d_{Lf}+d_{Rf}},
\qquad
\lambda_f^H
=
\left[
\frac{\theta_f}{\lambda_L}
+\frac{1-\theta_f}{\lambda_R}
\right]^{-1},
```

where $d_{Lf}$ and $d_{Rf}$ are the cell-center-to-face distances. Let $\mathcal{G}_f$ denote the positive orthogonal
area-over-distance factor in the over-relaxed face decomposition from Eq. {eq}`chap3_scalar_diffusion_flux_split_final`. The
two-point conductive energy rate, positive from $L$ to $R$, is

```{math}
F_{\lambda,f}^{\mathrm{orth}}
=
\mathcal{G}_f\lambda_f^H(T_L-T_R).
```

The same $\lambda_f^H$ weights the lagged nonorthogonal correction. The common-temperature closure assumes no unresolved thermal
contact resistance, so this Fourier connection remains active across a non-phase-changing material interface. On a link crossing
a reconstructed phase-change interface, it is replaced by the two one-sided heat fluxes in
Eq. {eq}`chap13_phase_change_energy_closure`.

## Governing Equations

### Species transport

For each owner material, the species conservation law advances
$z_{mk}=\alpha_m\rho_m\widetilde Y^{(m)}_k$ directly; $z_{mk}$ is the
segregated species unknown and $\widetilde Y^{(m)}_k$ is derived for a
positive material block:

```{math}
:label: chap13_species_transport_target

\frac{\partial z_{mk}}{\partial t}
+\nabla\cdot(z_{mk}\mathbf{u})
=-\nabla\cdot\boldsymbol{\mathcal{J}}^{(m)}_k
+\alpha_m\dot{\omega}^{(m)}_k
+S^{tr}_{mk}
```

for $k\in\mathcal{K}_m$. Here
$\dot{\omega}^{(m)}_k$ is an intrinsic production rate per unit volume of
material $m$, while $S^{tr}_{mk}$ is a conservative transfer or external source
per unit mixture volume. This distinction is important: intrinsic chemistry
vanishes with $\alpha_m$, but an inter-material transfer source must be able to
create mass in a receiver material that was absent before the transfer.

The mixture-area diffusion flux in Eq. {eq}`chap13_species_transport_target` is

```{math}
:label: chap13_mixture_area_diffusion_flux

\boldsymbol{\mathcal{J}}^{(m)}_k
=\alpha_m\mathbf{J}^{(m)}_k,
```

where $\mathbf{J}^{(m)}_k$ is the intrinsic diffusion flux per area occupied
by material $m$.

#### Species diffusion in material $m$

Molecular diffusion starts from a mixture-averaged Fick flux for the species
in material $m$. For each species owned by that material, form the uncorrected
intrinsic flux

```{math}
:label: chap13_material_raw_fick_flux

\widehat{\mathbf{J}}^{(m)}_k
=
-\rho_mD_{mk}\nabla\widetilde{Y}^{(m)}_k,
\qquad k\in\mathcal{K}_m,
```

where $D_{mk}$ is the mixture-averaged molecular diffusivity in
$\mathrm{m^2/s}$. Thus $\rho_mD_{mk}$ has units of
$\mathrm{kg/(m\,s)}$ and $\widehat{\mathbf{J}}^{(m)}_k$ has units of
$\mathrm{kg/(m^2\,s)}$. Independently evaluated mixture-averaged Fick fluxes
do not generally sum to zero. The material species equation nevertheless
requires

```{math}
:label: chap13_diffusion_constraint_internal

\sum_{k\in\mathcal{K}_m}\mathbf{J}^{(m)}_k=\mathbf{0},
\qquad
\sum_{k\in\mathcal{K}_m}\boldsymbol{\mathcal{J}}^{(m)}_k=\mathbf{0}.
```

The selected CMIM method imposes this constraint on the integrated raw normal
flux at each finite-volume face by the sign-group rescaling defined below. It
does not additionally apply a composition-weighted correction velocity. The
correction is performed separately for each material, so it never balances one
material's flux against another material's flux. A one-species material has
identically zero corrected species diffusion. Inter-material mass transfer is
supplied by $S^{tr}_{mk}$, not by passive Fick diffusion.

For mass-conserving chemistry within material $m$,

```{math}
:label: chap13_intrinsic_chemistry_mass_constraint

\sum_{k\in\mathcal{K}_m}\dot{\omega}^{(m)}_k=0.
```

The transfer sources satisfy

```{math}
:label: chap13_material_transfer_source_constraint

\sum_{k\in\mathcal{K}_m}S^{tr}_{mk}=\dot{m}_m,
\qquad
\sum_{m=1}^{N_m}\sum_{k\in\mathcal{K}_m}S^{tr}_{mk}
=\dot{m}_{tot}.
```

### Material and mixture mass balances

Summing Eq. {eq}`chap13_species_transport_target` over the species of
material $m$ gives

```{math}
:label: chap13_material_continuity_from_species

\frac{\partial}{\partial t}
\sum_{k\in\mathcal K_m}z_{mk}
+\nabla\cdot\left(\mathbf u\sum_{k\in\mathcal K_m}z_{mk}\right)
=
-\nabla\cdot\sum_{k\in\mathcal K_m}\boldsymbol{\mathcal J}^{(m)}_k
+\alpha_m\sum_{k\in\mathcal K_m}\dot\omega_k^{(m)}
+\sum_{k\in\mathcal K_m}S_{mk}^{tr}.
```

The diffusion and chemistry sums vanish within each material. With
$q_m=\sum_k z_{mk}$ and $S_m\equiv\dot m_m=\sum_k S_{mk}^{tr}$, the result is

```{math}
:label: chap13_material_mass_transport

\frac{\partial q_m}{\partial t}
+\nabla\cdot(q_m\mathbf u)=S_m.
```

This balance follows from the species equations. Diffusion changes the species
distribution within a material, but contributes no material mass flux. The EOS
density then gives $\alpha_m=q_m/\rho_m$ for interface reconstruction.

Summing once more over materials, with $\rho=\sum_m q_m$ and
$\dot m_{tot}=\sum_m S_m$, gives mixture continuity:

```{math}
:label: chap13_continuity

\frac{\partial\rho}{\partial t}
+\nabla\cdot(\rho\mathbf u)=\dot m_{tot}.
```

Closed inter-material transfer conserves the integral of $\dot m_{tot}$.
For Kunkelmann--Hardt--Wondra redistribution, removal and addition occur in
different cells, so the local source remains in continuity even though its
integral over a closed domain is zero.

### Mixture momentum

```{math}
:label: chap13_momentum

\frac{\partial(\rho\mathbf{u})}{\partial t}
+\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})
=-\nabla p+\nabla\cdot\boldsymbol{\tau}
+\rho\mathbf{g}+\mathbf{f}_{\sigma}+\mathbf{S}_u
```

Here $\boldsymbol{\tau}$ is the viscous stress tensor, $\mathbf{g}$ is the
gravitational acceleration, $\mathbf{f}_{\sigma}$ is the surface-tension force
density, and $\mathbf{S}_u$ is any additional momentum source.


### Mixture static enthalpy

Define the material-resolved species-enthalpy diffusion flux

```{math}
:label: chap13_cmim_species_enthalpy_flux

\mathbf{q}_{hY}
=
\sum_{m=1}^{N_m}
\sum_{k\in\mathcal{K}_m}
\bar h_{mk}\boldsymbol{\mathcal{J}}^{(m)}_k
=
\sum_{m=1}^{N_m}
\alpha_m
\sum_{k\in\mathcal{K}_m}
\bar h_{mk}\mathbf{J}^{(m)}_k,
```

where $\bar h_{mk}$ is the partial specific enthalpy supplied by material $m$'s EOS
state. Both the enthalpy and diffusion flux retain their material index; this
notation does not imply a single global EOS or a global cross-material
diffusion flux.

The static-enthalpy balance is

```{math}
:label: chap13_enthalpy

\frac{\partial(\rho h)}{\partial t}
+\nabla\cdot(\rho\mathbf{u}h)
=\frac{Dp}{Dt}
+\nabla\cdot(\lambda\nabla T)
+\Phi
-\nabla\cdot\mathbf{q}_{hY}
+S_{\rho h}
```

Here $Dp/Dt$ is pressure work, $\nabla\cdot(\lambda\nabla T)$ is Fourier
conduction, $\Phi$ is viscous dissipation, and
$-\nabla\cdot\mathbf{q}_{hY}$ is enthalpy transport by the corrected physical
species-diffusion operator. The conservative static-enthalpy source
$S_{\rho h}$ is derived from the same transferred-energy terms as
$S_{\rho h_t}$ after kinetic-energy and mechanical-work bookkeeping. It is not
an additional chemistry or latent-heat source to be added independently. If
this conservative balance is rewritten for the specific variable $h$, the
corresponding $-h\dot{m}_{tot}$ correction must be retained.

For the convention in Eq. {eq}`chap13_momentum`, where $\mathbf{S}_u$ is the
additional conservative momentum source and gravity and capillarity are shown
separately, the coupled source terms satisfy

```{math}
:label: chap13_static_total_enthalpy_source_relation

S_{\rho h}
=
S_{\rho h_t}
-\mathbf{u}\cdot\mathbf{S}_u
+K\dot m_{tot}.
```

Thus transferred momentum, transferred kinetic energy, and thermal or chemical
energy cannot be specified independently. If a source contribution is moved
out of $\mathbf{S}_u$ or $S_{\rho h_t}$ and written explicitly, its paired term
must be moved at the same time.

Temperature is obtained from the common-temperature enthalpy inversion in
Eq. {eq}`chap13_temperature_from_enthalpy`, using the solved conservative
material variables $q_m^{*,k}$ and the material compositions derived from
$z_{ms}^{*,k}$. The material volume fractions are then derived from the EOS
densities at the corresponding closure state; they are not predictor unknowns.
The derived mixture species-fraction vector from
Eq. {eq}`chap13_global_from_internal` does not replace the independently
evaluated material thermodynamic states.

### Mixture total enthalpy

The mixture total enthalpy is defined by
Eq. {eq}`chap13_total_enthalpy_def`. Its conservative continuum balance is

```{math}
:label: chap13_cmim_total_enthalpy_balance

\frac{\partial(\rho h_t-p)}{\partial t}
+\nabla\cdot(\rho\mathbf{u}h_t)
=
\nabla\cdot(\lambda\nabla T)
-\nabla\cdot\mathbf{q}_{hY}
+\nabla\cdot(\boldsymbol{\tau}\cdot\mathbf{u})
+\rho\mathbf{g}\cdot\mathbf{u}
+\mathbf{f}_{\sigma}\cdot\mathbf{u}
+S_{\rho h_t}.
```

The conserved energy variable is $\rho h_t-p=\rho E$. The selected
$h_t$ equation places the pressure-time contribution on the right-hand side;
that rearrangement does not introduce a new physical source. The conservative
source $S_{\rho h_t}$ contains any remaining energy
transfer, including energy carried by transferred mass. If continuity has a
mass source and the equation is rewritten for the specific variable $h_t$,
the corresponding $-h_t\dot{m}_{tot}$ product-rule correction must be retained.

The mechanical-work terms must be consistent with the forces retained in the
momentum equation. If a force is represented in momentum but its work is
omitted from total enthalpy, that omission is an additional modeling
approximation rather than part of the species-diffusion closure. In particular,
gravity and surface-tension work are retained whenever the corresponding
forces are retained in momentum.

For fixed control volumes with BDF1 time integration, the pressure-time
rearrangement has the following cell form:

```{math}
:label: chap13_total_enthalpy_bdf1_pressure_time

\frac{\Omega_P}{\Delta t}
\left[
(\rho h_t)_P^{n+1}-(\rho h_t)_P^n
\right]
+\sum_f F^{adv}_{\rho h_t,P,f}
=
\frac{\Omega_P}{\Delta t}
\left(p_P^{n+1}-p_P^n\right)
+\mathcal{R}_{E,P},
```

where $F^{adv}_{\rho h_t,P,f}$ is positive outward from $P$, and
$\mathcal{R}_{E,P}$ contains the consistently signed conductive,
diffusive-enthalpy, viscous-work, body-force-work, capillary-work, and
conservative source contributions. The pressure used by this term is lagged
during a segregated iteration, and outer convergence recovers the displayed
time-discrete balance.

On a fixed boundary with outward normal $\mathbf{n}$, the outward total-energy
flux associated with Eq. {eq}`chap13_cmim_total_enthalpy_balance` is

```{math}
:label: chap13_total_energy_boundary_flux

\mathbf{F}_E
=
\rho\mathbf{u}h_t
-\lambda\nabla T
+\mathbf{q}_{hY}
-\boldsymbol{\tau}\cdot\mathbf{u}.
```

Its advective contribution is evaluated from the material-resolved fluxes, not
from an unrelated mixture interpolation. An impermeable, adiabatic, no-work
wall sets the corresponding normal mass, heat, species-diffusion, and
mechanical-work fluxes to zero. Open boundaries must prescribe incoming
material fractions and material compositions together with a thermodynamically
compatible thermal state. Any imposed diffusive species flux carries the same
$\bar h_{mk}$ multiplier used in $\mathbf{q}_{hY}$.

#### Discrete enthalpy flux due to species diffusion

At an internal face, the energy equation reuses the species flux from
Eq. {eq}`chap13_material_corrected_face_flux`:

```{math}
:label: chap13_cmim_species_enthalpy_face_flux

F_{hY,f}^{k}
=
\sum_{m:\,\mathscr S_{m,f}}
\sum_{s\in\mathcal{K}_m}
\bar h_{ms,f}^{k}F^{k}_{ms,f}.
```

This is the same physical diffusion flux evaluated at the current outer
iterate and used in the species equation. An unsupported material link makes
zero contribution because its material-species diffusion flux is zero. Cell
species enthalpies are evaluated for active materials and are not gated by an
individual face's diffusion support.
On a supported link, the face enthalpy is the distance-centered thermodynamic
multiplier

```{math}
:label: chap13_cmim_species_enthalpy_face_value

\bar h_{ms,f}^{k}
=(1-\theta_f)\bar h_{ms,L}^{k}
+\theta_f\bar h_{ms,R}^{k},
\qquad s\in\mathcal K_m.
```

It is not selected by the advective upwind direction. Positive $F_{hY,f}$
leaves the left cell and enters the right cell, exactly as the underlying
species fluxes do.

Static species enthalpy is sufficient because the material-block flux is
zero-sum:

```{math}
:label: chap13_cmim_static_species_enthalpy_identity

\sum_{k\in\mathcal{K}_m}
\left(\bar h_{mk}+\frac{1}{2}\lVert\mathbf{u}\rVert^2\right)
\mathbf{J}^{(m)}_k
=
\sum_{k\in\mathcal{K}_m}
\bar h_{mk}\mathbf{J}^{(m)}_k.
```

For the same reason, adding one common enthalpy-reference offset to every
species in a material block does not change that material's diffusive energy
flux. This cancellation does not make the references of independent material
EOS models interchangeable for phase change, where mass transfers between
material blocks.

### Constitutive and term definitions

Viscous stress, $\boldsymbol{\tau}$, is defined as:

```{math}
:label: chap13_viscous_stress

\boldsymbol{\tau}
=
\mu\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)
-\frac{2}{3}\mu\left(\nabla\cdot\mathbf{u}\right)\mathbf{I}
```

Viscous dissipation in Eq. {eq}`chap13_enthalpy` is:

```{math}
:label: chap13_viscous_dissipation

\Phi=\boldsymbol{\tau}:\nabla\mathbf{u}
```

### Surface tension and CELESTE curvature

Surface tension is the tendency of an interface to reduce its area. For the
interface $\Gamma_{mn}$ between immiscible materials $m$ and $n$, the surface
tension coefficient $\sigma_{mn}$ is both a force per unit length and an
interfacial energy per unit area. Its units may therefore be written as either
$\mathrm{N/m}$ or $\mathrm{J/m^2}$. The coefficient does not depend on the order
of the pair, so $\sigma_{mn}=\sigma_{nm}$. The relations below use a constant
$\sigma_{mn}$. In this chapter, *capillary force* and *surface-tension force*
refer to the same force; capillarity describes the effects of surface tension
acting on a curved interface.

In VOF literature, a *color function* is a dimensionless scalar that identifies
the materials on the two sides of an interface. It is not a visible color, a
species fraction, or another thermodynamic state. Here it is obtained from the
material volume fractions. For the pair $m,n$, define the oriented pair color
function

```{math}
:label: chap13_celeste_pair_color

C_{mn}
=
\frac{\alpha_m}{\alpha_m+\alpha_n},
\qquad
\alpha_m+\alpha_n>0.
```

Thus $C_{mn}=1$ in material $m$ and $C_{mn}=0$ in material $n$. When these are
the only two materials, $\alpha_m+\alpha_n=1$ and $C_{mn}=\alpha_m$. If other
materials occupy part of a cell, the denominator makes $C_{mn}$ the relative
volume fraction within only the selected pair. Where
$\alpha_m+\alpha_n=0$, neither member of the pair is present and the local
volume fractions do not define a pair color.

Because the volume fractions are cell averages, a cell with
$0<C_{mn}<1$ contains both materials and is an interface cell. The nearby cells
over which the color changes from zero to one form the numerical interfacial
region, also called the interfacial band.

The word *oriented* states which material is assigned the value one. The
interface normal is a unit vector perpendicular to the local interface.
$\mathbf n_{mn}$ points in the direction of increasing $C_{mn}$, from material
$n$ toward material $m$. The curvature convention is

```{math}
:label: chap13_celeste_curvature_fit

\mathbf n_{mn}
=
\frac{\nabla C_{mn}}{\lVert\nabla C_{mn}\rVert},
\qquad
\kappa_{mn}=-\nabla\cdot\mathbf n_{mn},
\qquad
\lVert\nabla C_{mn}\rVert>0.
```

Curvature measures how strongly the interface bends and has units of inverse
length. With this orientation, a convex circular or spherical region of
material $m$ surrounded by material $n$ has positive curvature: the normal
points from the surrounding material into material $m$.

CELESTE (Curvature Evaluation with LEast-Squares fit of Taylor Expansion) [4]
evaluates the derivatives in Eq. {eq}`chap13_celeste_curvature_fit` from local
Taylor fits. In this context, *reconstruction* means estimating derivatives from
nearby cell values. CELESTE does not construct an explicit interface line,
plane, or surface. Let $P$ denote a cell, and let $\mathcal S_P$ be the
reconstruction stencil: the set of neighboring cells whose values are included
in the fit. Define

```{math}
:label: chap13_celeste_stencil_displacement

\mathbf s_{PQ}=\mathbf x_Q-\mathbf x_P,
\qquad Q\in\mathcal S_P.
```

Using $C$ as a short form of $C_{mn}$ in the local fit, the color gradient is
obtained from the second-order expansion

```{math}
:label: chap13_celeste_color_taylor_fit

C_Q-C_P
=
s_{PQ,i}\left.\frac{\partial C}{\partial x_i}\right|_P
+\frac{1}{2}s_{PQ,i}s_{PQ,j}
\left.\frac{\partial^2 C}{\partial x_i\partial x_j}\right|_P
+\mathcal O\!\left(\lVert\mathbf s_{PQ}\rVert^3\right).
```

Writing these equations as $\mathbf A_P\mathbf q_P\simeq\mathbf b_P$, where
$\mathbf q_P$ contains the gradient and the independent Hessian components,
the distance-weighted least-squares problem is

```{math}
:label: chap13_celeste_color_least_squares

\mathbf q_P
=
\underset{\mathbf q}{\operatorname{arg\,min}}
\left\lVert
\mathbf W_P\left(\mathbf A_P\mathbf q-\mathbf b_P\right)
\right\rVert_2^2,
\qquad
(\mathbf W_P)_{QQ}=\frac{1}{\lVert\mathbf s_{PQ}\rVert^2}.
```

Only the gradient components of $\mathbf q_P$ are required to form
$\mathbf n_{mn,P}$. For each component $a$ of this normal, a first-order Taylor
fit gives

```{math}
:label: chap13_celeste_normal_taylor_fit

n_{a,Q}-n_{a,P}
=
s_{PQ,j}G_{aj,P}
+\mathcal O\!\left(\lVert\mathbf s_{PQ}\rVert^2\right),
\qquad
G_{aj,P}=\left.\frac{\partial n_a}{\partial x_j}\right|_P.
```

The unweighted least-squares solution of Eq.
{eq}`chap13_celeste_normal_taylor_fit` supplies the preliminary curvature

```{math}
:label: chap13_celeste_raw_curvature

\kappa^{(0)}_P=-\operatorname{tr}(\mathbf G_P).
```

The raw curvature becomes unreliable as the color approaches zero or one,
where the color gradient vanishes. Two normalized averages extend a
representative curvature through the interfacial band. Define
$\mathcal S_P^+=\{P\}\cup\mathcal S_P$ and the color weight

```{math}
:label: chap13_celeste_color_weight

w_C(C)=\left[1-2\left|\frac{1}{2}-C\right|\right]^8,
\qquad 0\le C\le 1.
```

This weight is one at $C=1/2$ and zero in pure cells with $C=0$ or $C=1$.
It therefore gives the greatest weight to cells near the middle of the
numerical interface. The first average is

```{math}
:label: chap13_celeste_color_average

\kappa^{(C)}_P
=
\frac{
\displaystyle\sum_{Q\in\mathcal S_P^+}w_C(C_Q)\kappa^{(0)}_Q
}{
\displaystyle\sum_{Q\in\mathcal S_P^+}w_C(C_Q)
}.
```

For $Q\ne P$, let $\widehat{\mathbf s}_{PQ}=\mathbf s_{PQ}/
\lVert\mathbf s_{PQ}\rVert$ and define the directional weight

```{math}
:label: chap13_celeste_directional_weight

w_n(P,Q)=\left|\mathbf n_P\cdot\widehat{\mathbf s}_{PQ}\right|^8,
\qquad
w_n(P,P)=1.
```

The CELESTE curvature is the second normalized average

```{math}
:label: chap13_celeste_final_average

\kappa_P
=
\frac{
\displaystyle\sum_{Q\in\mathcal S_P^+}
w_C(C_Q)w_n(P,Q)\kappa^{(C)}_Q
}{
\displaystyle\sum_{Q\in\mathcal S_P^+}w_C(C_Q)w_n(P,Q)
}.
```

Both averages preserve a constant curvature. Reversing the pair orientation
gives

```{math}
:label: chap13_celeste_pair_orientation

\mathbf n_{nm}=-\mathbf n_{mn},
\qquad
\kappa_{nm}=-\kappa_{mn}.
```

The color gradient also changes sign when the pair orientation is reversed.
Consequently, the surface-tension force below is independent of which member
of the pair is named first.

At a sharp interface, surface tension acts on a surface rather than throughout
the surrounding volume. The continuum-surface-force (CSF) model [1] spreads
this surface force over the finite-width numerical interfacial region. Let
$\delta_{\Gamma,mn}$ be the resulting interfacial area per unit volume. It has
units of inverse length, and its volume integral gives the interface area,
$A_{\Gamma,mn}=\int_\Omega\delta_{\Gamma,mn}\,dV$. In the continuum-surface-
force representation,

```{math}
:label: chap13_surface_delta_color_gradient

\mathbf n_{mn}\delta_{\Gamma,mn}=\nabla C_{mn}.
```

The surface-tension force density is therefore

```{math}
:label: chap13_surface_tension_force

\mathbf f_{\sigma,mn}
=
\sigma_{mn}\kappa_{mn}\mathbf n_{mn}\delta_{\Gamma,mn}
=
\sigma_{mn}\kappa_{mn}\nabla C_{mn},
\qquad
\mathbf f_\sigma
=
\sum_{m<n}\mathbf f_{\sigma,mn}.
```

Here $\mathbf f_{\sigma,mn}$ is a force per unit volume. Each material pair
appears once in the sum. For a circular or spherical region of material $m$
with radius $R$, $C_{mn}=1$ inside the region and

```{math}
:label: chap13_surface_tension_laplace_jump

\kappa_{mn}=\frac{d-1}{R},
\qquad
p_m-p_n=\sigma_{mn}\kappa_{mn},
```

where $d=2$ for a circular cylinder and $d=3$ for a sphere. A *balanced-force*
discretization evaluates the pressure gradient and the capillary force with
compatible discrete operators. The pressure jump and surface-tension force can
then cancel in every cell at Young--Laplace equilibrium. A mismatch leaves a
small numerical force that drives nonphysical velocity near an otherwise static
interface. These velocities are called parasitic or spurious currents.

#### Interface contact with a solid wall

Where a material interface meets a solid wall, their intersection is the
*contact line*; it is a point in a two-dimensional section and a line in three
dimensions. The curved interface next to the wall is commonly called a
*meniscus*.

The *contact angle* is the angle between the wall and the tangent to the
interface, measured through a specified material. Naming that material is
essential: the angle measured through the other material is its supplement,
$180^\circ-\theta$. At equilibrium, the contact angle follows from the balance
of the three interfacial tensions at the contact line, known as Young's
relation.

A prescribed static contact angle fixes the interface direction at the wall and
therefore supplies a boundary condition for the interface normal and curvature.
An angle smaller than $90^\circ$ identifies the material as wetting the wall;
an angle larger than $90^\circ$ identifies it as non-wetting. In a narrow tube
or channel, the curved meniscus creates a Young--Laplace pressure difference.
The balance of this pressure difference with gravity produces capillary rise
for a wetting material or capillary depression for a non-wetting material.

A static contact angle specifies the equilibrium geometry but does not specify
how quickly the contact line moves. A dynamic contact-angle model may make the
angle depend on contact-line motion and may distinguish advancing and receding
angles. The contact-angle convention and the pair-color orientation must remain
consistent when the wall normal is used in the curvature calculation.

#### Surface-tension work in the total-enthalpy equation

The surface-tension force in the momentum equation does work on the fluid. The
conservative total-energy balance represented by the total-enthalpy equation
therefore uses the same force dotted with velocity. For a finite-volume cell
$P$, this contribution to the cell-integrated equation is

```{math}
:label: chap13_surface_tension_cell_work

\int_{\Omega_P}\mathbf f_\sigma\cdot\mathbf u\,dV
\approx
\mathbf u_P\cdot
\left(\int_{\Omega_P}\mathbf f_\sigma\,dV\right)
=
\Omega_P\mathbf f_{\sigma,P}\cdot\mathbf u_P.
```

Here $\mathbf f_{\sigma,P}$ is the cell-average force density. The force in
parentheses is the same cell-integrated force that enters the momentum
equation. It is not reconstructed independently for the energy equation. The
velocity in the work term is evaluated consistently with that momentum force,
including any time centering.

For a closed interface with constant $\sigma$ and no contact-line work,

```{math}
:label: chap13_surface_tension_work

\int_\Omega\mathbf{f}_\sigma\cdot\mathbf{u}\,dV
=
-\frac{d}{dt}
\left(\sum_{m<n}\sigma_{mn}A_{\Gamma,mn}\right).
```

Positive surface-tension work therefore increases the fluid total energy while
the interfacial energy decreases by the same amount. The sum of fluid total
energy and interfacial energy is conserved in a closed, otherwise source-free
system. Interfacial energy is determined by the interface area in this
constant-$\sigma$ model; it is not introduced as another transported bulk
state. Spatially varying surface tension or contact-line work requires the
corresponding additional energy terms.

### Source terms and sign conventions

Term definitions and sign conventions used in Eqs. {eq}`chap13_continuity`,
{eq}`chap13_material_mass_transport`, {eq}`chap13_species_transport_target`,
{eq}`chap13_momentum`, {eq}`chap13_enthalpy`, and
{eq}`chap13_cmim_total_enthalpy_balance`:

1. $\dot{m}_m$ is the net material mass source per mixture volume for material
   $m$; positive values add material mass.
2. $\dot{\omega}^{(m)}_k$ is an intrinsic species-production rate per material
   volume and is multiplied by $\alpha_m$.
3. $S^{tr}_{mk}$ is a conservative inter-material or external species source
   per mixture volume and is not multiplied by $\alpha_m$ again.
4. $S_{\rho h_t}$ contains conservative non-flux energy contributions,
   including energy carried by transferred mass and interfacial phase-change heat transfer.
   $S_{\rho h}$ is derived from it after kinetic and mechanical bookkeeping;
   chemistry or phase-change heat partitions must not be added to both.
5. $\mathbf{f}_{\sigma}$ is volumetric surface-tension force density.
6. $\dot{m}_{tot}$ is the local total mixture mass source. For closed surface
   phase change, its domain integral is zero but its local value generally is
   not.

## Species Discretization

The species equations determine both material transport and the composition
within each material. Volume-fraction reconstruction sets a common material
face flux, while material-local composition reconstruction distributes that
flux among the species.

### Face fluxes and cell balances

Following the solver framework in
{numref}`Chapter %s <chap_solution_algorithm>` and the finite-volume derivation in
{numref}`Chapter %s <chap_discretized_equations_unstructured>`, consider one
control volume $P$. The conserved species variable $z_{ms}$ is mass per
mixture volume. Its advective flux through a face is therefore a volume flow
rate times a species density:

```{math}
:label: chap13_species_advective_face_flux

F_{z,ms,f}=\dot{\Omega}_f z_{ms,f}
```

Here $\dot{\Omega}_f$ is the face volume flow rate defined in
Eq. {eq}`chap4_face_flux_definitions`, with units of volume per time.
Consequently $F_{z,ms,f}$ has units of mass per time. The generic scalar
equation in Eq. {eq}`chap3_discrete_scalar_transport_equation` instead uses
$\dot m_f\phi_f$, because its transported variable $\phi$ is a quantity per
unit mass. Both forms express advective transport through the face.

The face-only quantities use the stored
$\operatorname{cl}(f)$-to-$\operatorname{cr}(f)$ orientation. Summing the
species fluxes within material $m$ gives its mass flow rate:

```{math}
:label: chap13_q_face_flux

F_{q,m,f}
=q_{m,f}\dot{\Omega}_f,
\qquad
\dot{\Omega}_f=\vec{u}_f\mathbin{\cdot}\vec{S}_f
```

The effective face velocity $\vec u_f$ comes from momentum interpolation;
$\vec S_f$ is the oriented face-area vector.
The incidence sign $\sigma_{Pf}$ from {ref}`app_loci_face_orientation` converts
these stored fluxes to the outward orientation of cell $P$:

```{math}
:label: chap13_q_cell_outward_flux

\dot{\Omega}_{Pf}=\sigma_{Pf}\dot{\Omega}_f,
\qquad
F_{q,m,Pf}=\sigma_{Pf}F_{q,m,f}
```

A positive $\dot{\Omega}_f$ means flow from $\operatorname{cl}(f)$ to
$\operatorname{cr}(f)$. It is outflow when $P=\operatorname{cl}(f)$ and inflow
when $P=\operatorname{cr}(f)$. In contrast, $\dot{\Omega}_{Pf}>0$ always means outflow
from $P$, and $\dot{\Omega}_{Pf}<0$ always means inflow. This sign convention selects
the upwind state and supplies the signed contribution to the cell balance.
Integrating the species conservation law,
Eq. {eq}`chap13_species_transport_target`, over the fixed cell and applying
BDF1 gives

```{math}
:label: chap13_species_bdf1_balance

\frac{\Omega_P}{\Delta t}
\left(z_{ms,P}^{n+1}-z_{ms,P}^{n}\right)
+\sum_{f\in\mathcal F_P}
\left[\sigma_{Pf}F_{z,ms,f}^{n+1}+F_{ms,P,f}^{n+1}\right]
=\Omega_P S_{z,ms,P}^{n+1}
```

The diffusion flux $F_{ms,P,f}$ is the outward integral of
$\mathcal J_s^{(m)}$ through the face, and
$S_{z,ms}=\alpha_m\dot\omega_s^{(m)}+S_{ms}^{tr}$ is the species source per
mixture volume. Summing over the species of material $m$ cancels its net
species diffusion flux and gives the material balance:

```{math}
:label: chap13_q_bdf1_balance

\frac{\Omega_P}{\Delta t}
\left(q_{m,P}^{n+1}-q_{m,P}^{n}\right)
+
\sum_{f\in\mathcal F_P}F_{q,m,Pf}^{n+1}
=S_{m,P}^{n+1}\Omega_P
```

For one internal face, the two cell-outward contributions are

```{math}
:label: chap13_q_internal_face_conservation

F_{q,m,\operatorname{cl}(f)f}=+F_{q,m,f},
\qquad
F_{q,m,\operatorname{cr}(f)f}=-F_{q,m,f}
```

The shared flux therefore enters the two cell balances with opposite signs.
This material balance is obtained by summing the converged species rows.
The following sections define their face reconstructions and linearization.

### Volume-fraction face reconstruction

#### FOU and SOU volume-fraction candidates

For the face-reconstruction formulas in this subsection, write
$L_f=\operatorname{cl}(f)$ and $R_f=\operatorname{cr}(f)$. The stored
$\vec S_f$ points from $L_f$ to $R_f$, and
$\dot{\Omega}_f=\vec u_f\mathbin{\cdot}\vec S_f$. The donor $D_f$ is
$L_f$ when $\dot{\Omega}_f>0$ and $R_f$ when $\dot{\Omega}_f<0$; the other cell is the
acceptor. For compactness, the reconstruction equations below write these
cells with the subscripts $D$ and $A$.
This notation constructs one shared face state. Equation
{eq}`chap13_q_cell_outward_flux` supplies its sign in the balance
for cell $P$. First-order upwinding (FOU) uses

```{math}
\alpha_{m,f}^{FOU}=\alpha_{m,D}.
```

Second-order upwinding (SOU) applies the scalar Barth limiter independently to
each material volume fraction. With $\vec d_{Df}=\vec x_f-\vec x_D$, the
tentative component slope is

```{math}
r_{m,f}
=
\psi_{m,D}\nabla\alpha_{m,D}\mathbin{\cdot}\vec d_{Df}.
```

The vector reconstruction applies its nonnegativity guards component by
component:

```{math}
s_{m,f}^{(1)}
=
\begin{cases}
-\alpha_{m,D}, & \alpha_{m,D}+r_{m,f}<0,\\
r_{m,f}, & \text{otherwise},
\end{cases}
\qquad
\bar s_{m,f}
=
\begin{cases}
0, & \alpha_{m,D}-\dfrac32s_{m,f}^{(1)}<0,\\
s_{m,f}^{(1)}, & \text{otherwise}.
\end{cases}
```

The raw SOU candidate is

```{math}
\alpha_{m,f}^{SOU,raw}=\alpha_{m,D}+\bar s_{m,f}.
```

If $\alpha_{m,D}=0$, the guards set $\bar s_{m,f}=0$ for that component.
Other components retain their own SOU candidates. The coupled projection below
then restores a zero sum of the complete slope vector.

#### MIG volume-fraction candidate

During the coupled iteration, the EOS-derived fractions can sum to a value
slightly different from one. MIG first forms temporary fractions that describe
a complete partition of each cell in the reconstruction stencil:

```{math}
\alpha_\Sigma=\sum_m\alpha_m
\qquad
\beta_m=\frac{\alpha_m}{\alpha_\Sigma}
```

Here $\alpha_\Sigma$ must be finite and positive. Boundary fractions are
normalized in the same way. The usual discrete cell-gradient operator is then
applied to these $\beta_m$ values to obtain $\nabla\beta_m$. Thus the values
and gradients describe the same normalized field throughout the stencil.
They satisfy $\sum_m\beta_m=1$ and, to roundoff,
$\sum_m\nabla\beta_m=\vec0$, without changing the stored masses or
EOS-derived cell fractions. Section {ref}`chap13_mig_consistent_reconstruction`
derives this construction and explains its effect on trace-material transport.

MIG then constructs a normalized-variable face value using a virtual upstream
value bounded to the physical volume-fraction interval. This spatial NVD
normalization is separate from the material normalization above.
With $\vec r_{DA}=\vec x_A-\vec x_D$, define

```{math}
:label: chap13_vof_normalized_variable

\beta_{m,U}^{v}
=
\min\!\left(1,\max\!\left(0
\beta_{m,A}
-2\nabla\beta_{m,D}\cdot\vec r_{DA}
\right)\right)
\qquad
\widehat{\beta}_{m,D}
=
\frac{\beta_{m,D}-\beta_{m,U}^{v}}
     {\beta_{m,A}-\beta_{m,U}^{v}}
```

The bound and the normalization use the same upstream value. This prevents an
out-of-range extrapolation from sustaining a finite outgoing fraction as the
donor fraction vanishes. If the denominator is zero or the normalized donor
falls outside $[0,1]$, use the donor value directly. Constant states also use
that fallback.

The selected normalized face mapping is defined below. Its raw volume-fraction
candidate is

```{math}
:label: chap13_vof_denormalized_face_value

\alpha_{m,f}^{MIG,raw}
=
\beta_{m,U}^{v}
+\widehat{\beta}_{m,f}^{MIG}
\left(\beta_{m,A}-\beta_{m,U}^{v}\right)
```

The modified-IG (MIG) reconstruction is the normalized-variable-diagram (NVD)
mapping of Thakur, Wright, and Neal [1]. Define the IG branch

```{math}
:label: chap13_vof_ig_mapping

\mathcal{I}(z)
=
\begin{cases}
z, & z<0\ \text{or}\ z>1,\\
3z-2z^2, & 0\le z<\tfrac12,\\
1, & \tfrac12\le z\le1
\end{cases}
```

and form the face gradient by opposite-volume weighting,

```{math}
:label: chap13_vof_mig_alignment

\nabla\beta_{m,f}
=
\frac{
\Omega_A\nabla\beta_{m,D}+\Omega_D\nabla\beta_{m,A}
}{\Omega_D+\Omega_A}
\qquad
\gamma_f
=
\sqrt{
\left|
\frac{\nabla\beta_{m,f}\cdot\vec r_{DA}}
{\lVert\nabla\beta_{m,f}\rVert\lVert\vec r_{DA}\rVert}
\right|
}
\qquad 0\le\gamma_f\le1
```

For a vanishing gradient or degenerate direction, set
$\gamma_f=1$. MIG blends the compressive IG branch with the normalized
identity mapping,

```{math}
:label: chap13_vof_mig_mapping

\mathcal{N}_{MIG}(z)
=
\gamma_f\mathcal{I}(z)+(1-\gamma_f)z
\qquad
\widehat\beta_{m,f}^{MIG}
=
\mathcal N_{MIG}(\widehat\beta_{m,D})
```

Alignment increases the contribution of the compressive IG branch;
misalignment moves the candidate toward the normalized identity. A degenerate
normalization uses the donor value, which also preserves a constant state. If
$q_{m,D}=0$, the raw MIG candidate is
$\alpha_{m,f}^{MIG,raw}=\beta_{m,D}=0$.

#### Coupled volume-fraction projection

Both SOU and MIG first produce componentwise candidates. Their slope
balancing uses the donor value supplied to that reconstruction: the original
fraction for SOU, and the temporary normalized fraction for MIG. Denote
this base value by

```{math}
a_{m,D}^{SOU}=\alpha_{m,D}
\qquad
a_{m,D}^{MIG}=\beta_{m,D}
```

For $X\in\{SOU,MIG\}$, define

```{math}
:label: chap13_vof_deferred_face_value

s_{m,f}^{X}
=
\alpha_{m,f}^{X,raw}-a_{m,D}^{X}
\qquad
\Sigma_f^{X,-}=\sum_m\min(0,s_{m,f}^{X})
\qquad
\Sigma_f^{X,+}=\sum_m\max(0,s_{m,f}^{X})
\qquad
\Sigma_f^{X}=\Sigma_f^{X,-}+\Sigma_f^{X,+}
```

The sign group responsible for a nonzero total slope is reduced:

```{math}
:label: chap13_vof_sign_group_rescaling

\delta\alpha_{m,f}^{X}
=
\begin{cases}
\displaystyle
s_{m,f}^{X}
\dfrac{\Sigma_f^{X,-}-\Sigma_f^{X}}{\Sigma_f^{X,-}}
& \Sigma_f^{X}<0\ \text{and}\ s_{m,f}^{X}<0,\\[8pt]
\displaystyle
s_{m,f}^{X}
\dfrac{\Sigma_f^{X,+}-\Sigma_f^{X}}{\Sigma_f^{X,+}}
& \Sigma_f^{X}>0\ \text{and}\ s_{m,f}^{X}>0,\\[8pt]
s_{m,f}^{X}, & \text{otherwise}
\end{cases}
```

Each active denominator is nonzero, and each applied factor lies in $[0,1]$.
The final reconstructed volume fraction is

```{math}
:label: chap13_vof_face_simplex_correction

\alpha_{m,f}^{X}=a_{m,D}^{X}+\delta\alpha_{m,f}^{X}
\qquad
\sum_m\delta\alpha_{m,f}^{X}=0
```

The projection only scales a candidate slope toward zero. It therefore
preserves each slope's sign, any component bounds already satisfied by the
base value and raw candidate, and the sum of the base values. An exact-zero
candidate increment remains zero. Therefore,

```{math}
:label: chap13_vof_face_sum

\sum_m\alpha_{m,f}^{FOU}
=\sum_m\alpha_{m,f}^{SOU}
=\sum_m\alpha_{m,D}
\qquad
\sum_m\alpha_{m,f}^{MIG}=\sum_m\beta_{m,D}=1
```

If the base values and raw candidates lie in $[0,1]$ and the base values
sum to one, the reconstructed face vector belongs to the material-volume
simplex. MIG meets this unit-sum condition through its temporary fractions;
FOU and SOU retain any incoming donor sum error. Below, $\alpha_{m,f}$
denotes the selected FOU, SOU, or MIG result.

#### Conditional material-state reconstruction

Material density, composition, and species enthalpies are defined only
where the material is present. For a present donor, let
$\mathcal S_m(D_f)=1$ when the material is also present throughout the donor's
one-ring gradient stencil. A conditional material scalar $\varphi_m$ then uses

```{math}
:label: chap13_material_local_sou_face_value

\varphi_{m,f}^{SOU}
=
\begin{cases}
\varphi_{m,D_f}
+\psi_{\varphi_m,D_f}\nabla\varphi_{m,D_f}\mathbin{\cdot}
(\vec x_f-\vec x_{D_f}),
& \mathcal S_m(D_f)=1,\\[4pt]
\varphi_{m,D_f},
& q_{m,D_f}>0\ \text{and}\ \mathcal S_m(D_f)=0.
\end{cases}
```

Composition uses the same support condition. Its SOU increments are projected
to zero sum within each material block, as defined in
Eq. {eq}`chap13_material_species_face_increment`. When
$q_{m,D_f}=0$, the reconstructed material mass is zero and no conditional
material state is used.

An unqualified face value denotes the selected reconstruction. A superscript
$k$ means that the value was reconstructed from outer iterate $k$.

### Material face mass density

The donor $D_f^k$ is selected by the sign of $\dot{\Omega}_f^k$. Volume fraction and
intrinsic density give the reconstructed material face density:

```{math}
:label: chap13_q_linearized_face_flux

q_{m,f}^{\mathrm{rec},k}
=
\begin{cases}
\rho_{m,f}^k\alpha_{m,f}^k, & q_{m,D_f^k}^k>0,\\
0, & q_{m,D_f^k}^k=0.
\end{cases}
```

FOU gives $q_{m,f}^{\mathrm{rec},k}=q_{m,D_f^k}^k$. Every species of
material $m$ uses the same reconstructed material face density.

### Material composition at a face

For SOU, form the componentwise limited slope for every species
in material $m$,

```{math}
:label: chap13_material_species_raw_slope

d_{ms,f}^{k}
=
\psi_{Y_{ms},D_f^{k}}^{k}
\nabla\widetilde Y_{ms,D_f^{k}}^{(m),k}\mathbin{\cdot}
(\vec x_f-\vec x_{D_f^{k}}),
\qquad s\in\mathcal K_m.
```

The vector-MUSCL positivity guards are applied component by component in two
steps:

```{math}
:label: chap13_material_species_guarded_slope

g_{ms,f}^{(1),k}
=
\begin{cases}
-\widetilde Y_{ms,D_f^{k}}^{(m),k},
& \widetilde Y_{ms,D_f^{k}}^{(m),k}+d_{ms,f}^{k}<0,\\[3pt]
d_{ms,f}^{k}, & \text{otherwise},
\end{cases}
\qquad
g_{ms,f}^{k}
=
\begin{cases}
0,
& \widetilde Y_{ms,D_f^{k}}^{(m),k}
  -\dfrac{3}{2}g_{ms,f}^{(1),k}<0,\\[5pt]
g_{ms,f}^{(1),k}, & \text{otherwise}.
\end{cases}
```

Define the slope sums within this material block. The outer-iteration
superscript is suppressed in the next two equations.

```{math}
:label: chap13_material_species_slope_sums

G_{m,f}=\sum_{s\in\mathcal K_m}g_{ms,f},
\qquad
G_{m,f}^{-}=\sum_{s\in\mathcal K_m}\min(0,g_{ms,f}),
\qquad
G_{m,f}^{+}=\sum_{s\in\mathcal K_m}\max(0,g_{ms,f}).
```

Only the sign group responsible for a nonzero total slope is reduced:

```{math}
:label: chap13_material_species_rescaled_slope

\widehat g_{ms,f}
=
\begin{cases}
\displaystyle
g_{ms,f}\frac{G_{m,f}^{-}-G_{m,f}}{G_{m,f}^{-}},
&G_{m,f}<0\ \text{and}\ g_{ms,f}<0,\\[9pt]
\displaystyle
g_{ms,f}\frac{G_{m,f}^{+}-G_{m,f}}{G_{m,f}^{+}},
&G_{m,f}>0\ \text{and}\ g_{ms,f}>0,\\[9pt]
g_{ms,f},&\text{otherwise}.
\end{cases}
```

The active denominator in either rescaling branch is nonzero, and
$\sum_{s\in\mathcal K_m}\widehat g_{ms,f}=0$. The selected Picard-lagged face
increment is therefore

```{math}
:label: chap13_material_species_face_increment

\delta\widetilde Y_{ms,f}^{k}
=
\begin{cases}
\widehat g_{ms,f}^{k},
& \text{conditional-state scheme is SOU and }
  \mathcal S_m(D_f^{k})=1,\\[4pt]
0, & \text{conditional-state scheme is FOU, or SOU is unsupported}.
\end{cases}
\qquad s\in\mathcal K_m.
```

The species-slope projection is applied independently to each material and
preserves that material's composition sum.

(chap13_separate_material_composition_corrections)=
### Species convection and deferred correction

Consider an internal face and write $D=D_f^k$ for the donor selected by the
incoming volume flow $\dot{\Omega}_f^k$. The following reconstruction uses
a material present in that donor, $q_{m,D}^k>0$. Its face composition is

```{math}
:label: chap13_species_correction_composition_increment

\widehat Y_{ms,f}^k
=\widetilde Y_{ms,D}^{(m),k}+\delta\widetilde Y_{ms,f}^k,
\qquad
\sum_{s\in\mathcal K_m}\delta\widetilde Y_{ms,f}^k=0
```

The increment is given by
Eq. {eq}`chap13_material_species_face_increment`. Multiplying this
composition by the reconstructed material density gives the reconstructed
species density and its advective flux:

```{math}
:label: chap13_reconstructed_species_flux

z_{ms,f}^{\mathrm{rec},k}
=q_{m,f}^{\mathrm{rec},k}\widehat Y_{ms,f}^k,
\qquad
F_{z,ms,f}^{\mathrm{rec},k}
=\dot{\Omega}_f^k z_{ms,f}^{\mathrm{rec},k}
```

#### From the reconstructed flux to the standard predictor

Add and subtract the incoming donor value inside the face flux:

```{math}
:label: chap13_species_donor_correction_identity

F_{z,ms,f}^{\mathrm{rec},k}
=\dot{\Omega}_f^k z_{ms,D}^k
+\dot{\Omega}_f^k
\left(z_{ms,f}^{\mathrm{rec},k}-z_{ms,D}^k\right)
```

This is an algebraic identity. The deferred-correction linearization makes
one further choice: replace the first donor value by the unknown
$z_{ms,D}^{*,k}$, while holding the bracket at the incoming iterate. The
standard predictor flux is then

```{math}
:label: chap13_standard_species_predictor_flux

F_{z,ms,f}^{\mathrm{std},*,k}
=\dot{\Omega}_f^k z_{ms,D}^{*,k}
+\dot{\Omega}_f^k
\left(z_{ms,f}^{\mathrm{rec},k}-z_{ms,D}^k\right)
```

The first term supplies implicit donor-cell convection. The second is the
*deferred correction*: a known mass flow during the linear solve. This is
the same split as Eq. {eq}`chap3_convection_split_sum`, applied to the
density variable $z_{ms}$ and its volume-flow coefficient. At convergence,
$z_{ms,D}^{*,k}=z_{ms,D}^k$, so the predictor flux equals the reconstructed
flux. With FOU reconstruction, the bracket is zero.

For cell $P$ and its neighbor $N_f$, the implicit part can be written using
the cell-outward volume flow:

```{math}
:label: chap13_standard_species_cell_flux

\begin{aligned}
\sigma_{Pf}F_{z,ms,f}^{\mathrm{std},*,k}
={}&\max\!\left(\dot{\Omega}_{Pf}^k,0\right)z_{ms,P}^{*,k}
+\min\!\left(\dot{\Omega}_{Pf}^k,0\right)z_{ms,N_f}^{*,k}\\
&+\sigma_{Pf}\dot{\Omega}_f^k
\left(z_{ms,f}^{\mathrm{rec},k}-z_{ms,D}^k\right)
\end{aligned}
```

In the species balance, Eq. {eq}`chap13_species_bdf1_balance`, an outgoing
flow therefore adds $\max(\dot{\Omega}_{Pf}^k,0)$ to the diagonal. Moving
an incoming neighbor term to the right gives the nonnegative neighbor
coefficient $-\min(\dot{\Omega}_{Pf}^k,0)$. The deferred term also moves
to the right, with the opposite sign. An outgoing correction reduces that
right-hand side; this is how a large deferred correction can produce a
negative species predictor even though donor-cell convection is implicit.

#### Material and composition parts of the correction

The deferred correction changes both the amount of material transported and
its species proportions. To separate these effects, substitute the face
composition and the identity
$z_{ms,D}^k=q_{m,D}^k\widetilde Y_{ms,D}^{(m),k}$:

```{math}
:label: chap13_species_correction_decomposition

\begin{aligned}
z_{ms,f}^{\mathrm{rec},k}-z_{ms,D}^k
&=q_{m,f}^{\mathrm{rec},k}
\left(\widetilde Y_{ms,D}^{(m),k}+\delta\widetilde Y_{ms,f}^k\right)
-q_{m,D}^k\widetilde Y_{ms,D}^{(m),k}\\
&=\left(q_{m,f}^{\mathrm{rec},k}-q_{m,D}^k\right)
\widetilde Y_{ms,D}^{(m),k}
+q_{m,f}^{\mathrm{rec},k}\delta\widetilde Y_{ms,f}^k
\end{aligned}
```

The first contribution changes material transport at the donor composition.
Its species sum is the material mass-flow correction:

```{math}
:label: chap13_species_correction_material_sum

\sum_{s\in\mathcal K_m}
\dot{\Omega}_f^k\left(q_{m,f}^{\mathrm{rec},k}-q_{m,D}^k\right)
\widetilde Y_{ms,D}^{(m),k}
=\dot{\Omega}_f^k\left(q_{m,f}^{\mathrm{rec},k}-q_{m,D}^k\right)
```

The second contribution changes the species shares of that flow. Its sum
within the material is zero:

```{math}
:label: chap13_species_correction_composition_sum

\sum_{s\in\mathcal K_m}
\dot{\Omega}_f^k q_{m,f}^{\mathrm{rec},k}\delta\widetilde Y_{ms,f}^k=0
```

Denote this composition contribution by $\delta F_{z,ms,f}^k$. Including
the absent-donor branch, its definition is

```{math}
:label: chap13_direct_species_face_correction

\delta F_{z,ms,f}^k
=\begin{cases}
\dot{\Omega}_f^k q_{m,f}^{\mathrm{rec},k}\delta\widetilde Y_{ms,f}^k,
& q_{m,D}^k>0\\
0, & q_{m,D}^k=0
\end{cases}
```

Thus $\delta F_z$ denotes only the composition correction, with

```{math}
:label: chap13_qz_face_correction_identity

\sum_{s\in\mathcal K_m}\delta F_{z,ms,f}^k=0
```

For a present donor, the standard flux can now be written as

```{math}
:label: chap13_standard_species_flux_decomposed

F_{z,ms,f}^{\mathrm{std},*,k}
=\dot{\Omega}_f^k z_{ms,D}^{*,k}
+\dot{\Omega}_f^k\left(q_{m,f}^{\mathrm{rec},k}-q_{m,D}^k\right)
\widetilde Y_{ms,D}^{(m),k}
+\delta F_{z,ms,f}^k
```

This is the same deferred flux, with its two explicit contributions written
separately. For example, a donor material flow of $10\ \mathrm{kg/s}$ at
species fractions $(0.75,0.25)$ carries $(7.5,2.5)\ \mathrm{kg/s}$.
Changing the material flow to $8\ \mathrm{kg/s}$ at that composition gives
$(6,2)\ \mathrm{kg/s}$. Changing its composition to $(0.60,0.40)$ then gives
$(4.8,3.2)\ \mathrm{kg/s}$. The composition correction is
$(-1.2,+1.2)\ \mathrm{kg/s}$ and leaves the material flow unchanged. Each
species flux enters the two neighboring cells with opposite signs.

Summing the standard species fluxes gives the material predictor flux:

```{math}
:label: chap13_standard_material_predictor_flux

F_{q,m,f}^{\mathrm{std},*,k}
=\sum_{s\in\mathcal K_m}F_{z,ms,f}^{\mathrm{std},*,k}
=\dot{\Omega}_f^k q_{m,D}^{*,k}
+\dot{\Omega}_f^k\left(q_{m,f}^{\mathrm{rec},k}-q_{m,D}^k\right)
```

Here $q_{m,D}^{*,k}=\sum_s z_{ms,D}^{*,k}$; there is no separate material
solve. At convergence, this gives
$F_{q,m,f}=\dot{\Omega}_f q_{m,f}^{\mathrm{rec}}$.

At an absent donor, both the incoming donor and reconstructed material
densities are zero. The entire deferred correction is set to zero without
evaluating a conditional composition. The implicit term
$\dot{\Omega}_f^k z_{ms,D}^{*,k}$ remains: material supplied to that donor
through another face can still pass through the linear stencil.

#### Optional changes to the convection split

Implicit material advection moves the material correction into the donor
coefficient, while retaining the explicit composition correction. Its
derivation and predictor flux are given in
{ref}`chap13_material_advection_variants`. The diffusion and cell-balance
construction below apply to either convection split.

Conservative rescaling, developed in {ref}`chap13_conservative_rescaling`,
can reduce either correction. Use $\widehat q_{m,f}^k$ for the resulting
material face density; without rescaling it equals
$q_{m,f}^{\mathrm{rec},k}$. The same selected density is used for every
species of material $m$. Rescaling also replaces $\delta F_{z,ms,f}^k$ by
the bounded composition correction in
Eq. {eq}`chap13_separate_species_face_correction`. The flux and pressure
relations below use these selected quantities.

### Finite-volume species diffusion

The generic finite-volume diffusion split is developed in
{ref}`chap3_scalar_transport_rows`. The closure below specializes its
compact implicit part and lagged nonorthogonal correction to the material
support and block conservation required by algebraic VOF.

The finite-volume operator combines intrinsic molecular transport with
algebraic-VOF material support in the cell coefficient

```{math}
:label: chap13_material_diffusion_coefficient

\Gamma_{mk}=\alpha_m\rho_mD_{mk}.
```

A material face link is supported only when

```{math}
:label: chap13_material_diffusion_support

\mathscr S_{m,f}
\iff
q_{m,L}>0\ \text{and}\ q_{m,R}>0.
```

Consider an internal face $f$ between cells $L$ and $R$, with positive flux
oriented from $L$ to $R$. Let $d_{Lf}$ and $d_{Rf}$ be the
cell-center-to-face distances and define

```{math}
:label: chap13_material_diffusion_face_weights

\theta_f=\frac{d_{Lf}}{d_{Lf}+d_{Rf}},
\qquad
1-\theta_f=\frac{d_{Rf}}{d_{Lf}+d_{Rf}}.
```

For unit face normal $\mathbf{n}_f$ oriented from $L$ to $R$, face area $A_f$,
and cell-center vector
$\mathbf d_f=\mathbf{x}_R-\mathbf{x}_L$, define the face area vector and
positive orthogonal geometry factor by

```{math}
:label: chap13_material_diffusion_geometry_factor

\mathbf S_f=A_f\mathbf n_f,
\qquad
\mathcal{G}_f
=
\frac{A_f}
{\mathbf{n}_f\cdot\mathbf d_f}.
```

The corresponding over-relaxed decomposition is

```{math}
:label: chap13_material_diffusion_geometry_split

\mathbf S_f
=
\mathcal G_f\mathbf d_f+\mathbf T_f,
\qquad
\mathbf T_f=\mathbf S_f-\mathcal G_f\mathbf d_f.
```

The distance-weighted harmonic face coefficient of the complete supported
coefficient is

```{math}
:label: chap13_material_diffusion_face_coefficient

\Gamma_{mk,f}^{H}
=
\begin{cases}
\displaystyle
\left[
\frac{\theta_f}{\Gamma_{mk,L}}
+\frac{1-\theta_f}{\Gamma_{mk,R}}
\right]^{-1},
& \mathscr S_{m,f},\ \Gamma_{mk,L}>0,\ \text{and}\ \Gamma_{mk,R}>0,\\[8pt]
0, & \text{otherwise}.
\end{cases}
```

The positive compact conductance is therefore

```{math}
:label: chap13_material_diffusion_face_conductance

G_{mk,f}=\mathcal G_f\Gamma_{mk,f}^{H}.
```

For a material containing only one species, the composition is identically
one and its numerical conductance and diffusion flux are therefore taken as
zero. Its sole direct $z_{mk}$ transport equation is still solved.

No small positive regularization is introduced. An unsupported link contributes
zero without using either conditional composition or its gradient in the
physical flux.

The ordinary conditional-composition gradient also has a support requirement
on its use. Let
$\mathscr R_{m,P}$ mean that material $m$ is present in cell $P$ and in every
sample used by the selected cell-gradient stencil. The lagged face gradient is
used when

```{math}
:label: chap13_material_diffusion_gradient_support

\mathscr R_{m,f}\iff \mathscr R_{m,L}\ \text{and}\ \mathscr R_{m,R},
```

in which case the shared opposite-cell-volume face weighting gives

```{math}
:label: chap13_material_diffusion_face_gradient

\overline{\nabla\widetilde Y}^{(m)}_{k,f}
=
\frac{
\Omega_R\nabla\widetilde Y^{(m)}_{k,L}
+\Omega_L\nabla\widetilde Y^{(m)}_{k,R}
}{\Omega_L+\Omega_R}.
```

Here $\Omega_L$ and $\Omega_R$ are the adjacent cell volumes.

The same material support is used for every species in the block. Consequently,
a normalized material composition has
$\sum_{k\in\mathcal K_m}\overline{\nabla\widetilde Y}^{(m)}_{k,f}=\mathbf0$
in exact arithmetic. When $\mathscr R_{m,f}$ is false but the two face-adjacent
cells still support the material, only the nonorthogonal correction is
suppressed; the compact two-point flux remains available.

For a supported link with a usable face gradient, substitute
Eq. {eq}`chap13_material_diffusion_geometry_split` into the integrated Fick
flux:

```{math}
\begin{aligned}
-\Gamma_{mk,f}^{H}
\nabla\widetilde Y^{(m)}_{k,f}\cdot\mathbf S_f
&\approx
-\Gamma_{mk,f}^{H}
\left[
\mathcal G_f
\left(\widetilde Y^{(m)}_{k,R}-\widetilde Y^{(m)}_{k,L}\right)
+\overline{\nabla\widetilde Y}^{(m)}_{k,f}\cdot\mathbf T_f
\right]\\
&=
G_{mk,f}
\left(\widetilde Y^{(m)}_{k,L}-\widetilde Y^{(m)}_{k,R}\right)
-\Gamma_{mk,f}^{H}
\overline{\nabla\widetilde Y}^{(m)}_{k,f}\cdot\mathbf T_f.
\end{aligned}
```

The complete raw integrated face flux is

```{math}
:label: chap13_material_raw_face_flux

\widehat F_{mk,f}
=
\begin{cases}
G_{mk,f}
\left(\widetilde Y^{(m)}_{k,L}-\widetilde Y^{(m)}_{k,R}\right)
-\Gamma_{mk,f}^{H}
\overline{\nabla\widetilde Y}^{(m)}_{k,f}\cdot\mathbf T_f,
& G_{mk,f}>0\ \text{and}\ \mathscr R_{m,f},\\[3pt]
G_{mk,f}
\left(\widetilde Y^{(m)}_{k,L}-\widetilde Y^{(m)}_{k,R}\right),
& G_{mk,f}>0\ \text{and}\ \neg\mathscr R_{m,f},\\
0, & G_{mk,f}=0.
\end{cases}
```

It approximates the flux of
$\boldsymbol{\mathcal{J}}^{(m)}_k$ through the face and has units of
$\mathrm{kg/s}$. It is positive from $L$ to $R$: the leading term is positive
when $\widetilde Y^{(m)}_{k,L}>\widetilde Y^{(m)}_{k,R}$. A diffusive source
defined as positive into cell $L$ has the opposite sign.

For one material block, define the positive and negative raw flux magnitudes

```{math}

P_{m,f}=\sum_{k\in\mathcal K_m}[\widehat F_{mk,f}]^+,
\qquad
N_{m,f}=\sum_{k\in\mathcal K_m}[-\widehat F_{mk,f}]^+,
```

where $[a]^+=\max(a,0)$. To enforce a zero material-block sum without changing
any raw-flux sign, define the sign-group factors

```{math}

s^+_{m,f}
=
\begin{cases}
N_{m,f}/P_{m,f},
& P_{m,f}>N_{m,f},\\
1, & P_{m,f}\le N_{m,f},
\end{cases}
\qquad
s^-_{m,f}
=
\begin{cases}
P_{m,f}/N_{m,f},
& N_{m,f}>P_{m,f},\\
1, & N_{m,f}\le P_{m,f}.
\end{cases}
```

The selected denominator is strictly positive whenever either quotient is
evaluated, so this definition needs no small flux regularization.
The corrected flux is

```{math}
:label: chap13_material_corrected_face_flux

F_{mk,f}
=
\begin{cases}
s^+_{m,f}\widehat F_{mk,f},
& \text{supported link and }\widehat F_{mk,f}>0,\\
s^-_{m,f}\widehat F_{mk,f},
& \text{supported link and }\widehat F_{mk,f}\le 0,\\
0, & \text{unsupported link}.
\end{cases}
```

The sign-group correction is applied once to the complete raw flux in
Eq. {eq}`chap13_material_raw_face_flux`, not separately to its compact and
nonorthogonal parts. Only the sign group with the larger total magnitude is
reduced. No raw-flux sign changes, no magnitude increases, and an exactly zero
raw species flux remains zero. In exact arithmetic the correction enforces

```{math}
:label: chap13_material_species_diffusion_flux_sum

\sum_{k\in\mathcal{K}_m}F_{mk,f}=0
```

on every face, with only floating-point summation roundoff in the
implementation. If a block has only one species, its corrected
flux is zero. The face contributes $-F_{mk,f}$ to the left-cell right-hand side
and $+F_{mk,f}$ to the right-cell right-hand side, so conservation does not
require a face-composition interpolation or any change to the transported cell
composition.

In the segregated species solve, only the current compact conductance is
treated implicitly. The incoming-Picard $q_m^k$ is frozen during every scalar
$z_{ms}$ solve. For a supported link, define the left- and right-side scalar
coefficients

```{math}
:label: chap13_material_diffusion_z_coefficients

d_{ms,L}^{k}=\frac{G_{ms,f}^{k}}{q_{m,L}^{k}},
\qquad
d_{ms,R}^{k}=\frac{G_{ms,f}^{k}}{q_{m,R}^{k}}.
```

The conductance has units of $\mathrm{kg/s}$, while q and z have units of
$\mathrm{kg/m^3}$. Consequently $d_{ms,L}^{k}$ and $d_{ms,R}^{k}$ have the
$\mathrm{m^3/s}$ units required by the direct-z matrix. No division is
performed on an unsupported link because its conductance is exactly zero.

The complete corrected physical flux $F_{ms,f}^{k}$ is evaluated from the
current derived composition, including its supported nonorthogonal correction
and material-block conservation correction. The compact q-normalized value at
the current state is removed through the source, and the complete physical
flux is added. The diffusion flux represented by the linearized species
equation is therefore

```{math}
:label: chap13_material_linearized_diffusion_flux

F_{ms,f}^{lin,*,k}
=
F_{ms,f}^{k}
+G_{ms,f}^{k}
\left[
\left(\frac{z_{ms,L}^{*,k}}{q_{m,L}^{k}}
-\frac{z_{ms,R}^{*,k}}{q_{m,R}^{k}}\right)
-
\left(\frac{z_{ms,L}^{k}}{q_{m,L}^{k}}
-\frac{z_{ms,R}^{k}}{q_{m,R}^{k}}\right)
\right],
\qquad s\in\mathcal K_m.
```

Equivalently, the left row receives diagonal $+d_{ms,L}^{k}$ and right-neighbor
coefficient $-d_{ms,R}^{k}$, while the right row receives diagonal
$+d_{ms,R}^{k}$ and left-neighbor coefficient $-d_{ms,L}^{k}$. With

```{math}
:label: chap13_material_diffusion_source_replacement

C_{ms,f}^{k}
=G_{ms,f}^{k}
\left(\frac{z_{ms,L}^{k}}{q_{m,L}^{k}}
-\frac{z_{ms,R}^{k}}{q_{m,R}^{k}}\right)-F_{ms,f}^{k},
```

the left and right right-hand sides receive $+C_{ms,f}^{k}$ and
$-C_{ms,f}^{k}$, respectively. Evaluating the assembled rows at the incoming
state therefore leaves $+F_{ms,f}^{k}$ outward from the left cell and
$-F_{ms,f}^{k}$ outward from the right cell.

The value in Eq. {eq}`chap13_material_linearized_diffusion_flux` is positive
from $L$ to $R$. Its cell-outward form is

```{math}
:label: chap13_material_diffusion_cell_outward_flux

F_{ms,P,f}^{lin,*,k}
=
\begin{cases}
+F_{ms,f}^{lin,*,k}, & P=L,\\
-F_{ms,f}^{lin,*,k}, & P=R.
\end{cases}
```

The same equal-and-opposite mapping applies to the converged physical flux
$F_{mk,f}$. Thus the diffusion term in a cell row is positive outward from that
cell, while the two adjacent rows receive opposite contributions.

The identity $q_m=\sum_s z_{ms}$ holds at every iterate. During the linear
solve, however, $z_{ms}^{*,k}/q_m^k$ uses a lagged denominator and is not the
predicted composition $z_{ms}^{*,k}/q_m^{*,k}$. As the outer iteration
converges, these ratios agree and
$F_{ms,f}^{lin,*,k}=F_{ms,f}^{k}$. The physical flux $F_{ms,f}^{k}$ is used in
the species-diffusion energy term. It remains block-conservative by
Eq. {eq}`chap13_material_species_diffusion_flux_sum`; the temporary implicit correction is
not a second constitutive flux.

Harmonic interpolation of $\alpha_m\rho_mD_{mk}$ defines the
algebraic-VOF support model. A face joining two different pure materials has
zero species conductance, while species may diffuse where the same material
has positive support in both adjacent cells. Scalar volume fractions do not
contain enough topology to determine whether two portions of the same material
inside neighboring mixed cells are geometrically connected. Coefficient
interpolation alone cannot remove that limitation.

#### Scope of molecular diffusion

The molecular-diffusion model uses mass-fraction-gradient,
mixture-averaged Fick diffusion. An impermeable wall has zero normal molecular
species conductance and flux. At an open boundary, prescribed composition
sets the advected inflow state; a nonzero molecular species flux requires a
separate boundary-flux closure and cannot be inferred from that inflow value
alone. A captured material interface is handled by the internal VOF-support
operator, while phase-changing transfer is an explicit conservative source
and not passive cross-material Fick diffusion. A phase-changing wall is
therefore not an impermeable wall and requires a coupled boundary mass,
species, and energy sources.

On a nonorthogonal mesh, Eq. {eq}`chap13_material_raw_face_flux` adds the
lagged correction only where the conditional composition gradient has complete
material support. A support-crossing or incomplete gradient stencil falls back
to the supported two-point flux. This avoids using a dormant composition as
continuation data, but it does not reconstruct the subcell material aperture or
an interface-normal boundary condition. Because overlapping positive VOF tails
in neighboring cells can activate a material link, the resulting diffusion
remains sensitive to the resolved interface thickness.

Mass-fraction-gradient Fick diffusion is not automatically predictive for
dense liquids or near-critical cryogenic mixtures. Applications involving
liquid oxygen, liquid hydrogen, liquid methane, or related mixtures use diffusivities valid for the relevant
phase, pressure, temperature, and composition and are restricted to regimes
in which the mass-fraction-gradient approximation is appropriate.


### Species predictor row

Use the standard convective flux in
Eq. {eq}`chap13_standard_species_predictor_flux`, or the optional implicit
material flux in Eq. {eq}`chap13_direct_species_face_flux`, together with
the linearized diffusion flux in
Eq. {eq}`chap13_material_linearized_diffusion_flux`. The predictor form of
the BDF1 species balance, Eq. {eq}`chap13_species_bdf1_balance`, is

```{math}
:label: chap13_discrete_direct_species_equation

{}^{(z)}\!R_{ms,P}^{*,k}
=
\frac{\Omega_P}{\Delta t}
\left(z_{ms,P}^{*,k}-z_{ms,P}^{n}\right)
+\sum_{f\in\mathcal F_P}
\left[
\sigma_{Pf}F_{z,ms,f}^{*,k}
+F_{ms,P,f}^{lin,*,k}
\right]
-\Omega_PS_{z,ms,P}^{k}
=0.
```

Here $S_{z,ms}=\alpha_m\dot\omega_s^{(m)}+S_{ms}^{tr}$, and
$F_{ms,P,f}^{lin,*,k}$ is positive outward from cell $P$. The equation is
written for every $s\in\mathcal K_m$, including the sole species of a
one-species material.

Collecting the implicit terms gives the usual cell equation:

```{math}
:label: chap13_species_algebraic_row

{}^{(z)}\!a_{ms,P}^{k}z_{ms,P}^{*,k}
=
\sum_{N\in\mathcal N(P)}
{}^{(z)}\!a_{ms,P,N}^{k}z_{ms,N}^{*,k}
+{}^{(z)}\!b_{ms,P}^{k}.
```

The diagonal contains the time term, outgoing implicit convection, and compact
diffusion. The neighbor coefficients contain incoming convection and compact
diffusion. Accepted histories, physical sources, explicit diffusion
corrections, and the negative cell-outward deferred convective corrections
enter the right-hand side. For standard advection these include both the
material and composition contributions; for implicit material advection only
the composition correction remains explicit.
All coefficients and explicit terms are evaluated from the same incoming
iterate and held fixed during the solve.

For a species under-relaxation factor $0<\alpha_z\le1$, divide the diagonal
by $\alpha_z$ and add
$(1-\alpha_z)\,{}^{(z)}\!a_{ms,P}^{k}z_{ms,P}^{k}/\alpha_z$
to the right-hand side, using the unrelaxed diagonal in this addition.
This changes the iteration, but not the converged equation.

#### Summed material balance and convergence

The material predictor is defined by
$q_{m,P}^{*,k}=\sum_s z_{ms,P}^{*,k}$. Summing the species balances gives

```{math}
:label: chap13_q_picard_balance

\frac{\Omega_P}{\Delta t}
\left(q_{m,P}^{*,k}-q_{m,P}^{n}\right)
+\sum_{f\in\mathcal F_P}\sigma_{Pf}F_{q,m,f}^{*,k}
+\sum_{s\in\mathcal K_m}\sum_{f\in\mathcal F_P}F_{ms,P,f}^{lin,*,k}
=\Omega_PS_{m,P}^{k}.
```

This is the sum of the unrelaxed species equations; a relaxed solve also
contains the corresponding relaxation terms until outer convergence. The
convective sum is exact at every iteration. The corrected physical diffusion
flux also sums to zero within each material, but its lagged implicit correction
need not do so before convergence. That correction vanishes when
$z_{ms}^{*,k}=z_{ms}^{k}$, recovering the material balance in
Eq. {eq}`chap13_q_bdf1_balance`.

Thus $q_m=\sum_s z_{ms}$ is an identity at every stage, whereas convergence
of the material balance requires convergence of the species iteration.
A small residual of the frozen linear system alone is not sufficient.
Along with the equation residuals, measure the conservative species update:

```{math}
:label: chap13_species_update

\max_P
\frac{\displaystyle\sum_m\sum_{s\in\mathcal K_m}
\left|z_{ms,P}^{*,k}-z_{ms,P}^{k}\right|}
{\max\!\left(\rho_P^{*,k},\rho_P^{k}\right)}
\le\tau_z.
```

Here $\rho_P^{*,k}=\sum_m q_{m,P}^{*,k}$ and $\tau_z$ is the species-update
tolerance. Scaling by mixture density avoids division by a vanishing species
or material mass. Thermodynamic admissibility and material-volume closure
must also hold.

## Pressure-Correction Equation

After the species predictor, the pressure-correction solve holds
$z_{ms}^{*,k}$ and hence $q_m^{*,k}=\sum_s z_{ms}^{*,k}$ fixed. The
material reconstruction coefficients and explicit composition corrections
also remain fixed. The pre-pressure common-temperature
closure solves Eq. {eq}`chap13_temperature_from_enthalpy` at $p_P^k$ while
$q_{m,P}^{*,k}$, $z_{ms,P}^{*,k}$, $h_{t,P}^{*,k}$, and
$\vec{u}_P^{*,k}$ remain fixed. Denote this root by
$T_P^{\mathrm{pred},k}$. For every active material,
$q_m=\alpha_m\rho_m$ gives $\alpha_m=q_m/\rho_m$. At the pre-pressure state,
this gives

```{math}
\rho_{m,P}^{\mathrm{pred},k}
=
\rho_m\!\left(
p_P^k,
T_P^{\mathrm{pred},k},
\widetilde{\vec{Y}}_{P}^{(m),*,k}
\right),
\qquad
\alpha_{m,P}^{\mathrm{pred},k}
=
\frac{q_{m,P}^{*,k}}{\rho_{m,P}^{\mathrm{pred},k}}.
```

For an absent material, $q_{m,P}^{*,k}=0$ and
$\alpha_{m,P}^{\mathrm{pred},k}=0$; no material density is evaluated.
The heat capacities and density derivatives carrying a $\mathrm{pred},k$
label below are evaluated at the same state. Pressure correction then acts on
the volumes occupied by the fixed material masses and on the common carrier
volume flux. The cell and face linearizations below follow those two
dependencies.

### Predicted material-volume closure error

The predicted volume fractions must fill the cell. Their closure error is

```{math}
:label: chap13_fixed_q_volume_error

E_{V,P}^{\mathrm{pred},k}
=
\sum_{\substack{m\\q_{m,P}^{*,k}>0}}
\alpha_{m,P}^{\mathrm{pred},k}
-1.
```

- $E_{V,P}^{\mathrm{pred},k}=0$: the predicted material volumes exactly fill
  the cell.

- $E_{V,P}^{\mathrm{pred},k}>0$: the predicted materials require more volume
  than the cell contains.

- $E_{V,P}^{\mathrm{pred},k}<0$: the predicted materials leave part of the
  cell volume unfilled.

This is the local error in $\sum_m\alpha_m=1$, not the continuity-equation
residual. Pressure correction changes the material densities, and therefore
their implied volumes, while q remains fixed.

Absent materials contribute no volume and require no intrinsic density.
The summed species predictor already contains accepted history, face
transport, and physical material sources, so
Eq. {eq}`chap13_fixed_q_volume_error` is the only cell material-volume closure
term on the pressure right-hand side.

For each active material, let $\rho_{m,p}$ be the EOS density derivative at
fixed temperature and material composition, $\rho_{m,T}$ its derivative at
fixed pressure and composition, and $c_{p,m}=h_{m,T}$ the fixed-pressure heat
capacity. With the species densities fixed, their material sums and compositions are
also fixed. Define

```{math}
:label: chap13_fixed_q_temperature_pressure_response

A_P^k
=
\sum_{\substack{m\\q_{m,P}^{*,k}>0}}
q_{m,P}^{*,k}c_{p,m}^{\mathrm{pred},k},
\qquad
B_P^k
=
\sum_{\substack{m\\q_{m,P}^{*,k}>0}}
q_{m,P}^{*,k}
\left[
\frac{1}{\rho_{m,P}^{\mathrm{pred},k}}
+
\frac{T_P^{\mathrm{pred},k}\rho_{m,T}^{\mathrm{pred},k}}
{\left(\rho_{m,P}^{\mathrm{pred},k}\right)^2}
\right],
\qquad
\left(\frac{dT}{dp}\right)_{h,z}
=-
\frac{B_P^k}{A_P^k}.
```

The bracket in $B_P^k$ is
$h_{m,p}=1/\rho_m+T\rho_{m,T}/\rho_m^2$. Thus $A_P^k$ and $B_P^k$ are the
temperature and pressure derivatives of the fixed-q mixture-enthalpy
condition. The fixed-enthalpy volume compressibility is

```{math}
:label: chap13_fixed_q_volume_response

\chi_P^k
\equiv
-\left.
\frac{d}{dp}
\left(
\sum_{\substack{m\\q_{m,P}^{*,k}>0}}
\frac{q_{m,P}^{*,k}}{\rho_m}
\right)
\right|_{h,z}
=
\sum_{\substack{m\\q_{m,P}^{*,k}>0}}
\frac{q_{m,P}^{*,k}}
{\left(\rho_{m,P}^{\mathrm{pred},k}\right)^2}
\left[
\rho_{m,p}^{\mathrm{pred},k}
+
\rho_{m,T}^{\mathrm{pred},k}
\left(\frac{dT}{dp}\right)_{h,z}
\right].
```

This differentiates the common-temperature enthalpy condition. The
linearization assumes positive heat capacity, $A_P^k>0$, and positive
volume compressibility, $\chi_P^k>0$. For BDF1, the temporal diagonal
contribution is

```{math}
:label: chap13_fixed_q_temporal_pressure_coefficient

C_P^t
=
a_0\gamma_p\Omega_P\chi_P^{k},
```

where $a_0$ is the BDF1 new-time coefficient and $\gamma_p$ is Stream's
algorithmic compressible pressure-correction factor.

### Face-flux response

With implicit material advection, the material convective flux is the sum
of the species fluxes in
Eq. {eq}`chap13_separate_species_predictor_sum`. Its explicit composition
correction has zero material sum. Holding the donor, reconstructed coefficient,
and predicted species densities fixed gives

```{math}
:label: chap13_fixed_q_face_flux_derivative

\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f}
=
\begin{cases}
\displaystyle
\frac{\widehat q_{m,f}^k}{q_{m,D_f^k}^k}q_{m,D_f^k}^{*,k},
& q_{m,D_f^k}^k>0,\\[6pt]
q_{m,D_f^k}^{*,k}, & q_{m,D_f^k}^k=0.
\end{cases}
```

This response applies to the implicit material linearization. For standard
deferred material advection, use
Eq. {eq}`chap13_advection_variant_pressure_response`. The first branch above
includes the same reconstructed material coefficient as the species matrix.
The second retains the implicit FOU response of an initially absent donor.
With those coefficients held fixed, the linearized material
flux at a changed volume flux is

```{math}
:label: chap13_fixed_q_frozen_correction_flux

F_{q,m,f}(\dot{\Omega}_f)
=
F_{q,m,f}^{*,k}
+\left(\dot{\Omega}_f-\dot{\Omega}_f^k\right)
\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f}.
```

For the pressure row of cell $P$, convert each active material's mass-flux
response into a volume-flux response using its predicted cell density:

```{math}
:label: chap13_fixed_q_receiver_face_weight

W_{P,f}^k
=
\sum_{\substack{m\\q_{m,P}^{*,k}>0}}
\frac{1}{\rho_{m,P}^{\mathrm{pred},k}}
\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f}.
```

The active-material set is fixed during this linear solve. An absent material
in cell $P$ has no density with which to form a volume response. If pressure
correction activates a zero-flux face or reverses its donor, the next outer
iteration must account for transport into any newly receiving material.

For receiver $R$ and material $m$, let
$\mathcal F^{\mathrm{app}}_{R,m}$ contain such incoming faces whose corrected
donors contain material $m$ while the predicted receiver has none. A
donor-based estimate of the volume-fraction change is

```{math}
:label: chap13_pressure_appearance_update

\Delta\alpha_{\mathrm{app},R,m}^{k}
=
\frac{\Delta t}{\Omega_R}
\sum_{f\in\mathcal F^{\mathrm{app}}_{R,m}}
|\dot{\Omega}_f^{\mathrm c,k}|
\frac{q_{m,D_f}^{*,k}}{\rho_{m,D_f}^{\mathrm{pred},k}}.
```

Here $D_f$ is the donor selected by the corrected flux. The estimate uses only
present donor states. Summing per receiver and material and comparing with
the material-update tolerance detects an appearance change that still needs
another outer iteration.

For the pressure row of cell $P$, momentum interpolation supplies the
cell-outward face correction

```{math}
:label: chap13_fixed_q_face_volume_flux_correction

\dot{\Omega}_{Pf}^{\prime,k}
=d_f^{\mathrm{vol}}
\left(p_P^{\prime,k}-p_{N_f}^{\prime,k}\right),
\qquad
d_f^{\mathrm{vol}}
=
\frac{(1+\kappa_f)\,a_f^{p}}
{\rho_f}>0
```

Here $a_f^p$ is the face pressure-response coefficient before the SIMPLEC
relaxation factor $1+\kappa_f$ is applied, and $\rho_f$ is the mixture
face density used by momentum interpolation. Division by $\rho_f$ makes
$d_f^{\mathrm{vol}}$ a volume-flow response per unit pressure difference.
The fixed-q face coefficient for row $P$ is
$K_{P,f}=W_{P,f}^{k}d_f^{\mathrm{vol}}$. No independent pressure derivative of
$\widehat q_{m,f}^{k}$ or face intrinsic density is added.

### Algebraic pressure-correction row

This subsection writes each face flux in the outward orientation of cell $P$.
The species predictor uses $\dot{\Omega}_{Pf}^k$, while the momentum predictor supplies
$\dot{\Omega}_{Pf}^{*,k}$. The difference between these fluxes contributes to the pressure right-hand
side through the same material-volume weight used in the matrix:

```{math}
:label: chap13_fixed_q_pressure_rhs

B_P
=
a_0\Omega_PE_{V,P}^{\mathrm{pred},k}
-\sum_{f\in\mathcal F_P}
W_{P,f}^{k}
\left(\dot{\Omega}_{Pf}^{*,k}-\dot{\Omega}_{Pf}^{k}\right).
```

The internal-face row is

:::{admonition} CMIM Pressure-Correction Row

```{math}
:label: chap13_fixed_q_pressure_row

\left(
C_P^t+\sum_{f\in\mathcal F_P}K_{P,f}
\right)p_P^{\prime,k}
-\sum_{N\in\mathcal N(P)}
\left[
\sum_{f\in\mathcal F(P,N)}K_{P,f}
\right]p_N^{\prime,k}
=B_P.
```

The cell closure error uses the summed species predictor and its EOS state.
The flux difference, diagonal, and neighbor coefficients use the same
linearized face response. Physical mass sources already enter through the
species predictors; a second material-volume source would count their
contribution twice.
:::

Impermeable no-slip, slip, and symmetry boundaries have zero normal volume
flux. Periodic faces use the corresponding internal-face coupling. At an open
boundary, the flux and its pressure response follow the prescribed boundary
condition.


### Pressure, Velocity, and Volume-Flux Corrections

After solving for $p^{\prime,k}$, pressure and velocity are corrected as

```{math}
:label: chap13_cmim_pressure_under_relaxation

p_P^{\mathrm c,k}
=
\operatorname{clip}\!\left(
p_P^{k}+\alpha_p p_P^{\prime,k};
p_{\min},p_{\max}
\right),
```

where $\alpha_p$ is the pressure under-relaxation factor and
$\operatorname{clip}(x;p_{\min},p_{\max})$ confines $x$ to the configured
pressure interval.

```{math}
:label: chap13_cmim_velocity_pressure_update

(u_i)_P^{\mathrm c,k}
=
(u_i^{*,k})_P
-
\frac{
\displaystyle\sum_{f\in\mathcal F_P}
p_f^{\prime,k}(\vec S_{Pf})_i
}{{}^{(u)}\!a_P}.
```

The face volume flux is corrected by the momentum contribution:

```{math}
:label: chap13_cmim_volume_flux_pressure_update

\begin{aligned}
\dot{\Omega}_f^{\mathrm c,k}
&=
\dot{\Omega}_f^{*,k}+\dot{\Omega}_f^{\prime,k}
\\
&=
\dot{\Omega}_f^{*,k}
+d_f^{\mathrm{vol}}
\left(
p_{\operatorname{cl}(f)}^{\prime,k}
-p_{\operatorname{cr}(f)}^{\prime,k}
\right).
\end{aligned}
```

The same coefficient $d_f^{\mathrm{vol}}$ therefore appears in both the
pressure-correction row and the corrected common volume flux. Only this
momentum-induced contribution is applied to the common volume flux. The
material-density response remains in the pressure-correction equation, and
the outer iteration recovers density consistency.

The post-pressure continuity check evaluates the summed species flux at the
corrected volume flux, with the donor and reconstruction coefficients still
held at their predictor values. Its internal-face mass flux is

```{math}
:label: chap13_corrected_q_mass_flux

\dot m_f^{\mathrm{c},k}
=
\sum_m F_{q,m,f}^{*,k}
+
\left(\dot{\Omega}_f^{\mathrm c,k}-\dot{\Omega}_f^k\right)
\sum_m\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f}.
```

This is the mass flux used to assess continuity after pressure correction.
Momentum and total-enthalpy predictors instead use the incoming-state fluxes
defined below.

## Consistent Transport Requirement

The species, momentum, and total-enthalpy equations use the same incoming
material face reconstruction. Their linearizations contain different
predictor unknowns, but at nonlinear convergence they recover the common mass
and energy fluxes below. Each internal-face flux enters its adjacent cell
balances with opposite signs.

Let $L_f=\operatorname{cl}(f)$ and $R_f=\operatorname{cr}(f)$.
The stored face orientation and volume flux follow
Eq. {eq}`chap13_q_face_flux`. The donor is $L_f$ for $\dot{\Omega}_f>0$ and
$R_f$ for $\dot{\Omega}_f<0$.

Each material composition is a separate unit-sum vector. Reconstructing
species as one global vector would preserve only one total sum, not the sum
within every material. Instead, the block-local composition reconstruction
and support condition in Eq. {eq}`chap13_material_species_face_increment`
are applied separately to each material.

```{figure} media/cmim_face_state_construction.svg
:alt: Six finite-volume cells showing a notional helium-nitrogen boundary, an internal face, and upward material transport.
:name: chap13_cmim_face_state_construction
:align: center
:width: 100%

Material-resolved construction of an internal advective face flux for
$\dot{\Omega}_f>0$. The lower mixed cell is the donor and the upper cell is the
acceptor. The VOF scheme reconstructs each $\alpha_{m,f}$ at the marked face,
whereas the illustrated FOU/MIG path takes $\rho_{m,f}$,
$\widetilde Y^{(m)}_{k,f}$, and species enthalpies from the donor material state. The
SOU paths replace those donor values with the support-aware extrapolations in
Eqs. {eq}`chap13_material_local_sou_face_value` and
{eq}`chap13_material_species_face_increment`.
The curved helium--nitrogen boundary visualizes the cell-average occupancies;
it does not represent an explicitly reconstructed geometric interface.
```

### Material and species fluxes

The material mass flux associated with the reconstructed state is

```{math}
:label: chap13_material_flux

\dot m_{m,f}\equiv F_{q,m,f}
=\widehat q_{m,f}\dot{\Omega}_f.
```

Without conservative rescaling, $\widehat q_{m,f}=\rho_{m,f}\alpha_{m,f}$. With
conservative limiting it is the blend in
Eq. {eq}`chap13_separate_material_face_value`. The species fluxes are

```{math}
:label: chap13_species_flux

F_{mk,f}^{conv}
=\dot m_{m,f}\widetilde Y_{k,f}^{(m)},
\qquad k\in\mathcal K_m.
```

Here $\widetilde Y_{k,f}^{(m)}$ is the donor composition plus the selected
zero-sum increment. Without conservative rescaling it equals $\widehat Y_{mk,f}$;
with conservative limiting it is given by
Eq. {eq}`chap13_limited_face_composition`. Consequently,

```{math}
:label: chap13_species_face_flux_sum

\sum_{k\in\mathcal K_m}\widetilde Y_{k,f}^{(m)}=1,
\qquad
\sum_{k\in\mathcal K_m}F_{mk,f}^{conv}=\dot m_{m,f}.
```

If the material face mass flux is zero, all its advective species fluxes are
zero and no conditional face composition is needed. At a finite predictor
stage, Eq. {eq}`chap13_separate_species_predictor_sum` gives the
corresponding sum of the linearized species fluxes.

For FOU/FOU and FOU/MIG, intrinsic density, composition, and species
enthalpies come from the donor cell. For SOU/SOU and SOU/MIG, their
extrapolations require full material support on the donor's gradient stencil.
Unsupported conditional quantities retain the donor value. These are
extrapolations of cell states, not additional EOS evaluations at a face.

### Mixture mass and momentum fluxes

The mixture mass flux is the sum over materials:

```{math}
:label: chap13_mixture_material_flux_identity

\dot m_f=\sum_m\dot m_{m,f}.
```

At outer iteration $k$, the incoming material face densities give

```{math}
:label: chap13_incoming_picard_mass_flux

\rho_f^k=\sum_m\widehat q_{m,f}^k,
\qquad
\dot m_f^k=\rho_f^k\dot{\Omega}_f^k.
```

These are the face density and mass flux used by the momentum and
total-enthalpy predictors. They depend on the incoming state, not on a
completed species predictor.

The advective momentum flux is

```{math}
:label: chap13_material_momentum_flux_identity

\mathbf F_{\rho\mathbf u,f}^{adv}
=\dot m_f\mathbf u_f^{adv}
=\sum_m\dot m_{m,f}\mathbf u_f^{adv}.
```

The velocity $\mathbf u_f^{adv}$ is the selected FOU or SOU momentum face
state. It need not equal the pressure-coupled velocity
$\mathbf u_f$ that defines the normal volume flux.

### Total-enthalpy flux

Species enthalpy is transported by the same incoming species mass flux used
in its species equation. For a present donor, evaluating
Eq. {eq}`chap13_direct_species_face_flux` at the incoming state gives

```{math}
:label: chap13_incoming_direct_species_flux

F_{z,ms,f}^k
=
\dot{\Omega}_f^k\widehat q_{m,f}^k\widetilde Y_{ms,D_f^k}^{(m),k}
+\delta F_{z,ms,f}^k.
```

For an absent donor, this flux is zero without evaluating a composition.
Let $\bar h_{ms,f}^k$ be the FOU or support-aware SOU reconstruction of
the partial specific enthalpy already defined for species $s$ in material
$m$. The advected kinetic energy is
$K_f^k=\tfrac12\lVert\mathbf u_f^{adv,k}\rVert^2$. The total-enthalpy
flux is

```{math}
:label: chap13_incoming_direct_z_enthalpy_flux

F_{h,f}^{\mathrm{CMIM},k}
=
\sum_m\sum_{s\in\mathcal K_m}F_{z,ms,f}^k\bar h_{ms,f}^k
+\dot m_f^k K_f^k.
```

The zero-sum composition correction can carry enthalpy because the species
enthalpies need not be equal. Its contribution must therefore remain in this
energy flux even though it cancels from the material mass-flux sum.

Equivalently, define the material total enthalpy at the face by the same
effective species composition:

```{math}
:label: chap13_material_face_total_enthalpy

h_{t,m,f}
=\sum_{s\in\mathcal K_m}\widetilde Y_{s,f}^{(m)}\bar h_{ms,f}+K_f.
```

The total-enthalpy flux then has the material-sum form

```{math}
:label: chap13_material_enthalpy_flux_identity

F_{\rho h_t,f}^{adv}
=\sum_m\dot m_{m,f}h_{t,m,f}.
```

With FOU reconstruction, the material composition and species enthalpies are
their donor values, giving

```{math}
:label: chap13_donor_advective_kinetic_energy

h_{t,m,f}^{FOU}
=h_{m,D_f}+\frac12\lVert\mathbf u_{D_f}\rVert^2.
```

For SOU, the material face enthalpy follows from the reconstructed species
enthalpies and effective face composition above. An independently
extrapolated mixture or material enthalpy does not in general give the same
species-weighted flux.

For $\dot m_f\ne0$, the equivalent mixture face value is

```{math}
:label: chap13_consistent_mixture_enthalpy_face

h_{t,f}
=
\frac{\sum_m\dot m_{m,f}h_{t,m,f}}{\sum_m\dot m_{m,f}}.
```

The segregated total-enthalpy equation retains its implicit FOU term. Let
$F_{h,f}^{base,k}$ be the ordinary FOU or SOU enthalpy flux evaluated at the
incoming state. The material-resolved correction and linearized flux are

```{math}
:label: chap13_material_enthalpy_deferred_correction

\delta F_{h,f}^k
=F_{h,f}^{\mathrm{CMIM},k}-F_{h,f}^{base,k},
\qquad
F_{h,f}^{lin,*,k}
=F_{h,f}^{base,lin,*,k}+\delta F_{h,f}^k.
```

The correction enters adjacent cell balances with opposite signs and remains
fixed during the enthalpy solve. At nonlinear convergence,
$F_{h,f}^{base,lin,*,k}=F_{h,f}^{base,k}$, leaving the species-weighted
total-enthalpy flux. Molecular diffusion carries its separate enthalpy flux
from Eq. {eq}`chap13_cmim_species_enthalpy_face_flux`.

## Interface Heat Transfer and Phase Change

Phase change is an inter-material mass transfer. It is not passive species
diffusion across a material interface, and it is not obtained merely by
selecting a different EOS state. CMIM uses two distinct kinds of rate model:

| Model | Physical rate basis | Numerical representation |
| --- | --- | --- |
| Kunkelmann surface [7] | Local saturation equilibrium and one-sided heat transfer | Hardt--Wondra redistributed surface flux [8] |
| Schrage surface [10] | Molecular transfer kinetics coupled to the same heat balance | Hardt--Wondra redistributed surface flux [8] |
| Zwart volumetric [9] | Pressure departure and assumed bubble population | Collocated donor/receiver sources |

These models answer different physical questions. The surface model is
appropriate when the interface is resolved and heat supply limits the rate.
The volumetric model represents unresolved cavitation growth and collapse.
Both use the same signed donor-to-receiver convention and ultimately supply
the material, species, volume, and energy equations described below.

### Heat-limited transfer at a reconstructed interface

For transfer channel $c$, let $d(c)$ and $r(c)$ be the donor and receiver
materials, and let $\mathbf n_c$ point from donor to receiver across the
interface $\Gamma_c$. If the interface has no surface mass storage,
conservation of mass gives the signed Stefan flux

```{math}
:label: chap13_phase_change_interface_mass_jump

j_c''
=
\rho_{d,c}^{\Gamma}
(\mathbf u_{d,c}^{\Gamma}-\mathbf u_c^{\Gamma})\cdot\mathbf n_c
=
\rho_{r,c}^{\Gamma}
(\mathbf u_{r,c}^{\Gamma}-\mathbf u_c^{\Gamma})\cdot\mathbf n_c.
```

$j_c''>0$ transfers mass from donor to receiver; $j_c''<0$ reverses the
transfer through the same channel. For a liquid donor and vapor receiver,
these signs normally mean evaporation and condensation, respectively.

The local-equilibrium surface model evaluates the interface temperature from
a saturation relation,

```{math}
:label: chap13_kunkelmann_interface_state

T_c^{\Gamma}
=
T_{sat,c}\!\left(
p_c^{\Gamma},
\widetilde{\vec{Y}}_{d,c}^{\Gamma},
\widetilde{\vec{Y}}_{r,c}^{\Gamma}
\right).
```

For the present pure-component model, $p_c^\Gamma$ is the common local
pressure. One available saturation relation is the Antoine correlation,

```{math}
:label: chap13_antoine_saturation_relation

\log_{10}\!\left(\frac{p_{sat}}{p_0}\right)
=A-\frac{B}{T+C},
\qquad
T_{sat}(p)
=\frac{B}{A-\log_{10}(p/p_0)}-C.
```

The coefficients, pressure scale $p_0$, temperature convention, and validity
interval form one correlation and must be used consistently.

The alternative is the integrated Clausius–Clapeyron approximation. Assuming
ideal vapor, negligible liquid specific volume, and constant reference latent
heat $L_{ref}$ gives $d\ln p_{sat}/dT=L_{ref}/(R_vT^2)$. Integrating from a
specified saturation point $(T_{ref},p_{ref})$ gives

```{math}
:label: chap13_clausius_clapeyron_saturation

p_{sat}(T)=p_{ref}\exp\!\left[\frac{L_{ref}}{R_v}
\left(\frac{1}{T_{ref}}-\frac{1}{T}\right)\right],
\qquad R_v=\frac{R_u}{M_v}.
```

Here $R_v$ is the vapor species gas constant. This relation passes exactly
through the reference point; $L_{ref}$ sets its temperature sensitivity.
It can use the same evaluator and inverse as Antoine by setting
$p_0=p_{ref}$, $B=L_{ref}/(R_v\ln 10)$, $A=B/T_{ref}$, and $C=0$.
Both choices have a user-specified validity interval. Neither determines an
EOS phase envelope or resolves metastable EOS roots.

A multicomponent
interface would additionally require a partial-pressure, fugacity, or
chemical-potential closure. Likewise, using one interface pressure neglects
the saturation-temperature shifts associated with capillary pressure and
vapor recoil. Surface tension may still act in the momentum equation; that
does not by itself add a Kelvin correction to the saturation relation.

Let $a\in\mathcal C_c$ identify a conserved constituent transferred by the
channel, and let $k_d(a)$ and $k_r(a)$ locate it in the two material species
sets. If $\xi_{a,c}^\Gamma$ is its fraction of the transferred mass, the
specific enthalpy change between phases is

```{math}
:label: chap13_kunkelmann_phase_enthalpy

L_c^{\Gamma}
=
\sum_{a\in\mathcal C_c}\xi_{a,c}^{\Gamma}
\left(
\bar h_{r(c)k_r(a)}^{\Gamma}
-\bar h_{d(c)k_d(a)}^{\Gamma}
\right),
\qquad
\sum_a\xi_{a,c}^\Gamma=1.
```

The two partial specific enthalpies must use a compatible absolute reference
for each mapped constituent. With $\mathbf q=-\lambda\nabla T$ and no energy
stored at the interface, the Stefan heat balance is

```{math}
:label: chap13_kunkelmann_stefan_balance

j_c''L_c^{\Gamma}
=
\mathbf q_{d,c}^{\Gamma}\cdot\mathbf n_c
-\mathbf q_{r,c}^{\Gamma}\cdot\mathbf n_c
=
-\lambda_{d,c}(\nabla T)_{d,c}^{\Gamma}\cdot\mathbf n_c
+\lambda_{r,c}(\nabla T)_{r,c}^{\Gamma}\cdot\mathbf n_c.
```

The imbalance of heat arriving from the two phases supplies the enthalpy
change of the transferred mass. Consequently, latent heat is not a second
energy source to add after this balance has determined the rate.

Following Kunkelmann, the interface is the $C_c=1/2$ contour of an oriented
donor--receiver color field. A mixed-cell temperature is not generally a pure
phase temperature, so the normal gradients are sampled from pure material on
the two sides. If $T_{d,c}$ and $T_{r,c}$ are those samples and $d_{d,c}$ and
$d_{r,c}$ are their positive distances from the interface, then

```{math}
:label: chap13_kunkelmann_one_sided_gradients

(\nabla T)_{d,c}^{\Gamma}\cdot\mathbf n_c
\simeq \frac{T_c^{\Gamma}-T_{d,c}}{d_{d,c}},
\qquad
(\nabla T)_{r,c}^{\Gamma}\cdot\mathbf n_c
\simeq \frac{T_{r,c}-T_c^{\Gamma}}{d_{r,c}},
```

and the sampled heat-limited rate is

```{math}
:label: chap13_kunkelmann_sampled_stefan_flux

j_c''
\simeq
\frac{
\lambda_{d,c}(T_{d,c}-T_c^{\Gamma})/d_{d,c}
+\lambda_{r,c}(T_{r,c}-T_c^{\Gamma})/d_{r,c}
}{L_c^{\Gamma}}.
```

Both sides contribute to the rate and either side may favor evaporation or
condensation. When no direct pure-material neighbor exists, the one-sided
thermal information may be passed through a bounded mixed-cell band while
remaining tied to the reconstructed interface normal and distance.

The equilibrium closure assumes negligible kinetic resistance. Schrage replaces
that assumption with a molecular transfer law and solves for the interface
temperature, as developed next.

### Schrage kinetics coupled to the heat balance

Consider one pure-liquid donor and its pure-vapor receiver. The interface has
one temperature and no energy storage. Start with the same one-sided heat
supplies used in Eq. {eq}`chap13_kunkelmann_sampled_stefan_flux`. For either
side $s\in\{d,r\}$, write the heat entering the interface as

```{math}
:label: chap13_surface_thermal_coefficients

q_{s,c}''(T)=G_{s,c}(\Theta_{s,c}-T)
\qquad
G_{s,c}=\left\langle\frac{\lambda_{s,c}}{d_{s,c}}\right\rangle
\qquad
\Theta_{s,c}=\frac{\left\langle\lambda_{s,c}T_{s,c}/d_{s,c}\right\rangle}{G_{s,c}}
```

Here $G_{s,c}$ is a conductance per interface area, in W/(m² K), and
$\Theta_{s,c}$ is the effective pure-phase sample temperature. The brackets
are the normalized sampling and mixed-band propagation weights. A single
sample gives $G=\lambda/d$ and $\Theta=T_s$, recovering the preceding heat
balance. With several samples, averaging $G$ and $G\Theta$ preserves their
linear heat-supply relation. These coefficients are propagated before applying
the local interface temperature; propagating an already evaluated heat flux
would retain the temperature of the cell where it was first sampled.

For turbulent flow, the conductivity in these samples includes the same
eddy contribution as the bulk energy equation:

```{math}
:label: chap13_surface_effective_conductivity

\lambda_{s,eff}=\lambda_s+\frac{\mu_t c_p}{Pr_t}.
```

$\mu_t$ is the dynamic eddy viscosity, $Pr_t$ the turbulent Prandtl number,
and $c_p$ the solver's cell mixture heat capacity. Sampling uses cells that
meet the pure-material cutoff, so this heat capacity approaches the selected
phase value there. Molecular conductivity remains material-specific. The
same effective conductivity must be used when returning heat to the samples
below. Laminar flow has $\mu_t=0$. This is an eddy-diffusivity approximation
to normal interface heat transfer; it does not add an interface wall function
or a model for suppression of turbulence at the interface.

CMIM uses the base Menter SST 2003 equations for $k$ and $\omega$. Their
face mass flux is formed with CMIM's face density and the incoming volume
flux. For the optional buoyancy correction, a material mixture needs its
own thermal expansion coefficient. At fixed material mass fractions $w_m$,
$1/\rho=\sum_m w_m/\rho_m$; differentiation at fixed pressure gives

```{math}
:label: chap13_sst_mixture_expansion

\beta=-\frac{1}{\rho}\left(\frac{\partial\rho}{\partial T}\right)_{p,w}
=\sum_m\alpha_m\beta_m,
\qquad
\beta_m=-\frac{1}{\rho_m}\left(\frac{\partial\rho_m}{\partial T}\right)_{p,Y_m}.
```

Each present material supplies its EOS derivative; absent materials contribute
zero. The existing SST split is retained: positive buoyancy production enters
the $k$ source, while negative production is an implicit sink. This uses
thermal stratification, not the material-density jump across the interface,
to estimate buoyancy production. It is a RANS closure assumption that needs
validation for the tank, not an interface turbulence model. The current SST
path is limited to one species per material because material-local turbulent
species diffusion has not yet been added.

Schrage [10] supplies a kinetic expression for the same signed mass flux:

```{math}
:label: chap13_schrage_mass_flux

j_{S,c}''(T)=\frac{2\chi_c}{2-\chi_c}
\frac{1}{\sqrt{2\pi R_{v,c}}}
\left[
\frac{p_{sat,c}(T)}{\sqrt{T}}-
\frac{p_{v,c}}{\sqrt{T_{v,c}}}
\right]
\qquad
R_{v,c}=\frac{R_u}{M_{v,c}}
```

$\chi_c$ is the accommodation probability, with $0<\chi_c\le1$;
$M_{v,c}$ is the vapor molecular mass. The first term represents molecules
leaving the liquid and the second represents incident vapor molecules.
Positive flux means evaporation. In the pure-vapor approximation,
$p_{v,c}=p_c^\Gamma$. A noncondensable gas would require a vapor partial
pressure or a more general thermodynamic closure.

Two approximations are available for the incident-vapor temperature:
$T_{v,c}=\Theta_{r,c}$ uses the sampled vapor temperature, whereas
$T_{v,c}=T$ uses the common interface temperature. They differ when the vapor
sample and interface have different temperatures. The sampled value is a
mesh-based estimate, so its sensitivity to the sampling distance must be
examined. The accommodation probability is also a model input to qualify for
the chosen fluid and conditions.

Equating this kinetic rate to the Stefan heat-balance rate gives one scalar
equation for $T_c^\Gamma$:

```{math}
:label: chap13_schrage_temperature_residual

F_c(T)=G_{d,c}(\Theta_{d,c}-T)+G_{r,c}(\Theta_{r,c}-T)
-j_{S,c}''(T)L_c(p_c^\Gamma,T)=0
```

The latent enthalpy $L_c$ is the difference of the two material enthalpies in
Eq. {eq}`chap13_kunkelmann_phase_enthalpy`, evaluated at each trial temperature
and the common mechanical pressure. The saturation pressure belongs in the
kinetic law; it does not replace that pressure in either material EOS.
If Clausius–Clapeyron is selected, its input $L_{ref}$ controls the saturation
curve, while the energy balance still uses the EOS difference $L_c(p,T)$.
An exact constant-latent comparison requires material enthalpies satisfying
$h_r(p,T)-h_d(p,T)=L_{ref}$ throughout the comparison range. Equal constant
heat capacities with a fixed enthalpy offset provide an analytic test pair.
Replacing only the denominator by a prescribed constant would leave the
transported enthalpy sources inconsistent and is therefore not an option.
A bracketed scalar solve uses the intersection of the saturation correlation's
and material models' valid temperature ranges. It requires positive latent
enthalpy and checks both temperature resolution and the remaining heat-balance
residual. A missing bracket or invalid thermodynamic state is a model failure,
not a reason to clip the temperature to a range endpoint.

Once $T_c^\Gamma$ is known, the common heat-balance expression supplies
$j_c''=(q_{d,c}''+q_{r,c}'')/L_c^\Gamma$. Its agreement with
Eq. {eq}`chap13_schrage_mass_flux` is a verification check. The interface area,
redistribution, and conservative source equations below then apply to either
surface closure. No second kinetic or latent-heat source is added.

### Returning interface heat to the sampled cells

The two heat supplies must also appear in the energy equations of the cells
that supplied the temperatures. This remains necessary when the supplies
cancel: if $q_{d,c}''=-q_{r,c}''\ne0$, heat passes from one phase to the other
although $j_c''=0$. A source formed only from their sum would lose this
exchange.

For an interface patch in cell $P$, write the sampled heat supply from side
$s\in\{d,r\}$ as a sum over its pure-material sample cells $Q$:

```{math}
:label: chap13_phase_change_sample_heat

q_{s,c,P}''=\sum_Q g_{s,c,PQ}(T_Q-T_{c,P}^{\Gamma})
```

Here $g_{s,c,PQ}\ge0$ includes the thermal conductivity divided by sampling
distance and the normalized sampling weight; its units are
$\mathrm{W\,m^{-2}\,K^{-1}}$. With one direct sample it is simply
$\lambda_s/d_s$. When sampling passes through several mixed cells, the
normalized face-area weights along each path multiply, and contributions
from different paths add. Their conductance sum and temperature-weighted
sum give $G_{s,c,P}$ and $G_{s,c,P}\Theta_{s,c,P}$ used by the interface closure.

Each contribution in Eq. {eq}`chap13_phase_change_sample_heat` removes that
same amount of heat from its sample cell. Multiplying by the interface area
and summing the patches that use cell $Q$ gives the thermal energy source:

```{math}
:label: chap13_phase_change_heat_deposition

\Omega_Q\widehat S_{E,\Gamma,c,Q}^{pc}
=\sum_{P,s}A_{\Gamma,c,P}g_{s,c,PQ}
(T_{c,P}^{\Gamma}-T_Q)
```

Positive source heats the sample cell. The weights used for sampling are
therefore also used to return heat, in the reverse direction. Summing over
all sample cells gives

```{math}
:label: chap13_phase_change_deposited_heat_balance

\sum_Q\Omega_Q\widehat S_{E,\Gamma,c,Q}^{pc}
=-\sum_P A_{\Gamma,c,P}(q_{d,c,P}''+q_{r,c,P}'')
=-\sum_P A_{\Gamma,c,P}j_{c,P}''L_{c,P}^{\Gamma}
```

The separate side contributions can be nonzero even when this domain sum is
zero. Heat deposition uses no division by mass-transfer rate and is distinct
from the mass-source redistribution described below.

During a Picard iteration, interface temperature and conductance are held
fixed in the energy solve. Define
$K_Q=\sum_{P,s}A_{\Gamma,c,P}g_{s,c,PQ}$, in W/K. At fixed pressure,
composition, and velocity, the local temperature increment is approximated
by $\delta T_Q=\delta h_Q/c_{p,Q}$. The thermal source consequently contributes
$K_Q/c_{p,Q}$ to the energy diagonal and the matching
$(K_Q/c_{p,Q})h_Q^{(k)}$ to the right-hand side, in addition to the physical
source evaluated at iteration $k$. These matching terms cancel at
convergence. The interface closure is reevaluated at the next Picard
iteration.

### Geometry held fixed during a nonlinear solve

Reconstructing the contour from every algebraic-VOF iterate can abruptly
change the cut cells and one-sided sampling stencil during a segregated solve.
CMIM therefore reconstructs the surface geometry from the accepted
beginning-of-step volume fractions and holds it fixed during that step's
nonlinear iterations:

```{math}
:label: chap13_phase_change_accepted_geometry

C_c^{g}=C_c(\boldsymbol{\alpha}^{n}),
\qquad
\Gamma_c^{g}=\{\mathbf x:C_c^{g}(\mathbf x)=1/2\},
\qquad
\Gamma_c^{g,k}=\Gamma_c^{g}.
```

The current iterate still supplies pressure, temperature, conductivity,
interface thermodynamics, and $j_c''$. Thus only the discrete geometry is
lagged; the physical rate is neither clipped nor relaxed. The geometry is
rebuilt after a step is accepted. This first-order geometric treatment
requires the time step to resolve interface motion, and temporal refinement
should show convergence of the integrated heat and mass transfer.

### Hardt--Wondra conservative redistribution

The reconstructed surface rate first appears as a sharp volumetric source.
Using the surface delta distribution $\delta_{\Gamma_c}$,

```{math}
:label: chap13_phase_change_sharp_volume_source

s_{\Gamma,c}=j_c''\delta_{\Gamma_c},
\qquad
s_{\Gamma,c,P}\simeq
j_{c,P}''\frac{A_{\Gamma,c,P}}{\Omega_P},
\qquad
\dot M_c=\int_{\Gamma_c}j_c''\,dA.
```

Here $A_{\Gamma,c,P}$ is the reconstructed area in cell $P$. Direct use of
this thin, large-amplitude source is numerically difficult. Evaporation and
condensation can occur on different parts of the same interface. Split the
signed source into two nonnegative magnitudes before redistribution:

```{math}
:label: chap13_phase_change_directional_sources

s_{\Gamma,c}^{+}=\max(s_{\Gamma,c},0)
\qquad
s_{\Gamma,c}^{-}=\max(-s_{\Gamma,c},0)
\qquad
\dot M_c^{\pm}=\int_\Omega s_{\Gamma,c}^{\pm}\,dV
```

The original source is $s_{\Gamma,c}^{+}-s_{\Gamma,c}^{-}$ and the net rate is
$\dot M_c^{+}-\dot M_c^{-}$. The positive direction transfers donor to receiver;
the negative direction reverses that transfer. Both rates may be nonzero when
their difference is zero.

Apply the Hardt--Wondra smoothing operator [8] independently to each magnitude:

```{math}
:label: chap13_hardt_wondra_smoothing

\widetilde s_c^{\pm}
-\nabla\cdot\left(\ell_c^2\nabla\widetilde s_c^{\pm}\right)
=s_{\Gamma,c}^{\pm}
\qquad
\nabla\widetilde s_c^{\pm}\cdot\mathbf n_{\partial\Omega}=0
```

The numerical length $\ell_c$ controls how widely the source is spread; it is
not a physical interface thickness. The homogeneous Neumann condition
preserves each domain integral. Separating the directions before smoothing
prevents cancellation from erasing the amount transferred in either direction.

The smoothed fields are cropped to pure donor and receiver support. Let
$\chi_{d,c}$ and $\chi_{r,c}$ be the resulting zero-or-one indicators. They use
the accepted beginning-of-step fractions $\boldsymbol{\alpha}^{n}$ and remain
fixed during the nonlinear solve, like the interface geometry. This prevents
an intermediate volume-closure error from switching a cell's source on and off
as it crosses the pure-material cutoff. Rates and smoothing still update each
iteration.

For each direction, normalize the cropped field on each material side. Define
its support integral and deposition weight by

```{math}
:label: chap13_hardt_wondra_normalization

I_{m,c}^{\pm}=\int_\Omega\chi_{m,c}\widetilde s_c^{\pm}\,dV
\qquad
w_{m,c}^{\pm}=\frac{\chi_{m,c}\widetilde s_c^{\pm}}{I_{m,c}^{\pm}}
\qquad m\in\{d(c),r(c)\}
```

Each active weight has units of inverse volume and integrates to one. An
exactly zero directional rate has zero deposition weights and requires no
division. For a nonzero rate, missing pure support or a vanishing support
integral makes redistribution undefined. The finite-volume Helmholtz integral
is checked against each prescribed directional rate before normalization.

Multiplying these weights by the prescribed rates gives the material sources:

```{math}
:label: chap13_hardt_wondra_material_sources

S_{d(c),c}^{pc}=-\dot M_c^{+}w_{d,c}^{+}+\dot M_c^{-}w_{d,c}^{-}
\qquad
S_{r(c),c}^{pc}=+\dot M_c^{+}w_{r,c}^{+}-\dot M_c^{-}w_{r,c}^{-}
```

Their integrals are $-\dot M_c$ and $+\dot M_c$, with
$\dot M_c=\dot M_c^{+}-\dot M_c^{-}$. For evaporation alone this recovers the
usual normalization $N_{m,c}=\dot M_c/I_{m,c}$. Equal evaporation and condensation
retain spatial sources wherever their deposition weights differ, even though
the net transfer is zero. Nearby distributions can overlap and cancel locally;
the method retains the spatial smoothing approximation set by $\ell_c$.

Cropping keeps removal and addition out of mixed interface cells. Normalization
preserves the rates obtained from the heat balance. The support indicators
classify raw accepted volume fractions without clipping or renormalizing them.

For constituent-resolved transfer, the same mapping supplies equal and
opposite species masses:

```{math}
:label: chap13_phase_change_source_closure

\dot M_{a,c}
=
\int_{\Gamma_c}j_c''\xi_{a,c}^{\Gamma}\,dA,
\qquad
\int_\Omega S_{d(c)k_d(a),c}^{pc}\,dV=-\dot M_{a,c},
\qquad
\int_\Omega S_{r(c)k_r(a),c}^{pc}\,dV=+\dot M_{a,c},
\qquad
\sum_{k\in\mathcal K_m}S_{mk}^{pc}=S_m^{pc},
\qquad
\int_\Omega\sum_m S_m^{pc}\,dV=0.
```

All unmapped species receive zero source from that channel. A step that would
remove more donor mass than is available is inadmissible; changing the rate or
repairing the volume fractions would hide the inconsistency rather than solve
it.

### Coupling the surface transfer to the CMIM equations

Because surface removal and addition occur in different cells,
$\dot m_{tot}^{pc}=\sum_mS_m^{pc}$ need not vanish pointwise, even though its
domain integral is zero. It therefore remains in mixture continuity. The
corresponding material-volume source is

```{math}
:label: chap13_phase_change_material_volume_source

S_V^{pc}=\sum_m\frac{S_m^{pc}}{\rho_m}.
```

For constant phase densities, one channel integrates to
$\dot M_c(1/\rho_{r,c}-1/\rho_{d,c})$. Evaporation can therefore expand volume
while conserving mass. In the discrete pressure equation, this volume change
enters through the species sources, their predicted material sums, and the
resulting volume-closure error. It is not added again as a separate pressure
source. Under the common-velocity model, phase transfer has
no separate momentum source: $\mathbf S_u^{pc}=\mathbf0$.

Energy coupling follows the same Stefan balance. With
$K_c^\Gamma=\tfrac12\lVert\mathbf u_c^\Gamma\rVert^2$, define the donor and
receiver transferred total-enthalpy rates by

```{math}
:label: chap13_phase_change_energy_rates

\dot H_{t,m,c}^{\Gamma}
=
\int_{\Gamma_c}j_c''
\sum_{a\in\mathcal C_c}\xi_{a,c}^{\Gamma}
\left(\bar h_{m k_m(a)}^{\Gamma}+K_c^\Gamma\right)\,dA,
\qquad m\in\{d(c),r(c)\}.
```

For redistribution with both directions present, form separate versions
$\dot H_{t,m,c}^{\Gamma,+}$ and $\dot H_{t,m,c}^{\Gamma,-}$ of these integrals,
using $\max(j_c'',0)$ and $\max(-j_c'',0)$ as the mass-flux weights. Deposit them
with the same weights used for mass:

```{math}
:label: chap13_phase_change_directional_energy

\widehat S_{E,d,c}^{pc}
=-\dot H_{t,d,c}^{\Gamma,+}w_{d,c}^{+}
 +\dot H_{t,d,c}^{\Gamma,-}w_{d,c}^{-}
\qquad
\widehat S_{E,r,c}^{pc}
=+\dot H_{t,r,c}^{\Gamma,+}w_{r,c}^{+}
 -\dot H_{t,r,c}^{\Gamma,-}w_{r,c}^{-}
```

This construction preserves each directional energy integral without dividing
by net mass transfer. Interface enthalpy may differ between evaporating and
condensing regions. Within each direction, the deposition still uses an
integrated energy rate and the normalized mass distribution.

Here $k_m(a)$ is the mapped local species index for constituent $a$ in
material $m$. Also let
$\dot Q_c^\Gamma=\int_{\Gamma_c}(\mathbf q_{d,c}^\Gamma-\mathbf q_{r,c}^\Gamma)
\cdot\mathbf n_c\,dA$. The Stefan condition gives
$\dot H_{t,r,c}^\Gamma-\dot H_{t,d,c}^\Gamma=\dot Q_c^\Gamma$. The conservative
energy sources consequently satisfy

```{math}
:label: chap13_phase_change_energy_closure

\int_\Omega \widehat S_{E,d,c}^{pc}\,dV=-\dot H_{t,d,c}^{\Gamma},
\qquad
\int_\Omega \widehat S_{E,r,c}^{pc}\,dV=+\dot H_{t,r,c}^{\Gamma},
\qquad
\int_\Omega \widehat S_{E,\Gamma,c}^{pc}\,dV=-\dot Q_c^\Gamma.
```

These three terms sum to zero in a closed adiabatic domain. The interface term
represents the two one-sided conductive heat supplies deposited in their
sample cells by Eq. {eq}`chap13_phase_change_heat_deposition`. The ordinary Fourier
connection is therefore omitted only on mesh links that cross the
reconstructed interface; same-side conduction remains active. This counts
interfacial heat transfer exactly once and avoids adding a duplicate latent
heat source.

Let $\widehat S_E^{pc}$ be the sum of the three energy contributions. If
$\mathbf g=-\nabla\Phi_g$, the total-enthalpy and static-enthalpy equations use

```{math}
:label: chap13_phase_change_total_energy_source

S_{\rho h_t}^{pc}
=
\widehat S_E^{pc}-\Phi_g\dot m_{tot}^{pc},
\qquad
S_{\rho h}^{pc}
=
S_{\rho h_t}^{pc}+K\dot m_{tot}^{pc}.
```

The potential-energy term vanishes when gravity is absent. The second relation
is Eq. {eq}`chap13_static_total_enthalpy_source_relation` with
$\mathbf S_u^{pc}=\mathbf0$.

The source above belongs to the conservative total-enthalpy equation. With
`PDEContinuityAdded: yes` (the default), Stream subtracts $h_t$ times the
continuity equation from that equation. Because redistributed phase change can
add or remove mixture mass in an individual cell, the source must undergo the
same subtraction:

```{math}
:label: chap13_phase_change_assembled_energy_source

S_{h_t,\mathrm{assembled}}^{pc}
= S_{\rho h_t}^{pc} - h_t\dot m_{tot}^{pc}.
```

For example, adding mass at specific enthalpy $h_{t,\Gamma}$ gives the
conservative source $h_{t,\Gamma}\dot m_{tot}^{pc}$ and the assembled source
$(h_{t,\Gamma}-h_t)\dot m_{tot}^{pc}$. Thus mass entering at the local enthalpy
does not create an artificial temperature change. The subtraction vanishes
where the local mixture mass source is zero. Global mass conservation alone
does not make it vanish in each cell.

With `PDEContinuityAdded: no`, the conservative source is used directly. This
choice changes equation assembly, not the physical heat deposited or the
conservative source ledgers. Applying a common constant shift to both material
enthalpy references must leave the temperature and pressure histories unchanged;
this is checked by the closed-box evaporation and condensation tests.

### Pressure-driven volumetric transfer

A volumetric model does not reconstruct an interface or evaluate one-sided
heat fluxes. It computes a local rate $\dot m_c^{vol}$ directly from cell
state. For liquid donor and vapor receiver,

```{math}
:label: chap13_volumetric_phase_change_sources

S_{d(c),c}^{vol}=-\dot m_c^{vol},
\qquad
S_{r(c),c}^{vol}=+\dot m_c^{vol}.
```

These collocated sources sum to zero pointwise in mixture continuity, but
their material-volume contribution remains
$\dot m_c^{vol}(1/\rho_r-1/\rho_d)$.

The Zwart--Gerber--Belamri model [9] is based on a simplified
Rayleigh--Plesset bubble-growth relation. Let $\alpha_d$ and $\alpha_r$ be the
liquid and vapor volume fractions, $\alpha_{nuc}$ an assumed nuclei fraction,
$R_B$ a representative bubble radius, and $p_{sat}(T)$ the saturation
pressure. With positive rate from liquid to vapor,

```{math}
:label: chap13_zwart_volumetric_rate

\dot m_c^{vol}
=
\begin{cases}

C_{vap}
\dfrac{3\alpha_{nuc}\alpha_d\rho_r}{R_B}
\sqrt{\dfrac{2}{3}\dfrac{p_{sat}(T)-p}{\rho_d}},
& p<p_{sat}(T),
\\[1.1em]

-C_{cond}
\dfrac{3\alpha_r\rho_r}{R_B}
\sqrt{\dfrac{2}{3}\dfrac{p-p_{sat}(T)}{\rho_d}},
& p>p_{sat}(T),
\\[0.6em]

0,
& p=p_{sat}(T).

\end{cases}
```

$C_{vap}$ and $C_{cond}$ are empirical coefficients. The nuclei term controls
growth when pressure is below saturation, while the existing vapor fraction
controls collapse when pressure is above saturation. Merkle and
Schnerr--Sauer models use different time-scale or interfacial-area assumptions
but fill the same mathematical role.

Zwart is an inertial cavitation closure, not the heat-limited surface balance
in Eq. {eq}`chap13_kunkelmann_stefan_balance`. In the present collocated
form, no separate conservative mixture-energy source is added. The material
enthalpy partition changes while mixture total energy is held fixed, and the
EOS temperature inversion supplies the thermal response. A thermally limited
volumetric model would require an additional, consistent closure. The Zwart
and Kunkelmann rates should not simply be added because their physics may
overlap.

The present CMIM surface model is intentionally narrow: one mapped pure
constituent, one shared interface pressure, a common velocity, no interfacial
storage, and no kinetic resistance. Multicomponent equilibrium,
capillary-shifted saturation, and kinetic interface models are extensions to
this closure. These limits are kept explicit so that invalid thermodynamic or
material states are reported rather than hidden by clipping the state or
massaging the phase-transfer rate.

(chap13_segregated_nonlinear_algorithm)=
## Segregated Nonlinear Algorithm

Each outer iteration follows the pressure-based procedure in
{ref}`chap13_outer_iteration`. The accepted species, momentum, and energy
states supply the time histories. All predictor coefficients are evaluated
from the incoming iterate $k$ and held fixed during their respective linear
solves.

1. Derive $q_m^k=\sum_s z_{ms}^k$, $\rho^k=\sum_m q_m^k$, and each
   present material's composition. Evaluate the material EOS states at the
   common $(p^k,T^k)$ and derive $\alpha_m^k=q_m^k/\rho_m^k$.
2. Evaluate body-force, capillary, chemistry, and inter-material transfer
   terms. Use mutually consistent species, momentum, and total-energy
   sources. Material mass sources are the species-source sums. For surface
   phase change, hold the interface geometry in
   Eq. {eq}`chap13_phase_change_accepted_geometry` fixed during the nonlinear
   solve.
3. Reconstruct the material volume fractions and conditional face states.
   Assemble each species equation with the selected material-advection
   linearization, the explicit composition correction, and the physical and
   implicit diffusion terms. If conservative limiting is used, evaluate its
   factors from the same incoming state and accepted histories.
4. Solve the species, momentum, and total-enthalpy predictors. The momentum
   and energy equations use the incoming mass flux in
   Eq. {eq}`chap13_incoming_picard_mass_flux`; the enthalpy correction uses
   Eq. {eq}`chap13_incoming_direct_z_enthalpy_flux`. No predictor requires
   another predictor's newly solved value.
5. Derive $q_m^{*,k}=\sum_s z_{ms}^{*,k}$ and the predicted compositions.
   With the species predictors, $h_t^{*,k}$, $\mathbf u^{*,k}$, and $p^k$
   fixed, solve Eq. {eq}`chap13_temperature_from_enthalpy` for
   $T^{\mathrm{pred},k}$. Evaluate the predicted EOS densities, material
   volumes, and fixed-enthalpy volume response.
6. Use the predicted volume error and the difference between the momentum
   and species predictor volume fluxes to assemble
   Eq. {eq}`chap13_fixed_q_pressure_row`. Solve for $p^{\prime,k}$.
7. Correct pressure, velocity, and volume flux. At the corrected pressure,
   hold $z_{ms}^{*,k}$ and
   $h_t^{*,k}-\tfrac12\lVert\mathbf u^{\mathrm c,k}\rVert^2$ fixed and
   solve again for the common temperature. Reevaluate material densities and
   derive $\alpha_m^{k+1}=q_m^{*,k}/\rho_m^{k+1}$.
8. Check conservation-equation residuals, the species update in
   Eq. {eq}`chap13_species_update`, material-volume closure, enthalpy closure,
   and continuity using Eq. {eq}`chap13_corrected_q_mass_flux`. Account for
   pressure-induced material appearance using
   Eq. {eq}`chap13_pressure_appearance_update`. Accept the time level when
   the nonlinear tolerances and state-admissibility conditions are satisfied;
   otherwise repeat with the corrected state as iterate $k+1$.

The species sums define material mass throughout this procedure. The outer
iteration couples transport to pressure, temperature, and the derived volume
fractions.

## Limiting Balances and Reference States

Several limiting states follow directly from the governing balances and give
useful consistency conditions for any discretization of the framework.

A constant material or species enthalpy-reference shift leaves the physical
solution unchanged when it is applied consistently to the EOS enthalpy,
conservative energy state and history, face fluxes, boundary data, and sources.
The material identity $q_m=\sum_s z_{ms}$ follows from the species sum.
Discrete reference invariance additionally requires consistent species and
energy time terms, fluxes, and sources. Lagged implicit corrections recover
these identities at nonlinear convergence.

For a source-free, force-free, adiabatic stationary contact, let the material
states be evaluated at constants $p_0$ and $T_0$, with constant
material compositions. Then

```{math}
:label: chap13_stationary_contact_state

\mathbf{u}=\mathbf{0},
\qquad
p=p_0,
\qquad
T=T_0,
\qquad
\partial_t\alpha_m=0,
\qquad
q_m=\alpha_m\rho_m(p_0,T_0,\widetilde{\vec{Y}}^{(m)})
```

is an exact equilibrium even when intrinsic material densities differ. The
captured contact must not generate velocity, pressure, or temperature
departures merely because $\alpha_m$ changes across the interface.

On a periodic or unbounded domain, the corresponding uniform-velocity contact
has the exact translating solution

```{math}
:label: chap13_translating_contact_state

\mathbf{u}=\mathbf{u}_0,
\qquad
p=p_0,
\qquad
T=T_0,
\qquad
\alpha_m(\mathbf{x},t)
=
\alpha_m^0(\mathbf{x}-\mathbf{u}_0t).
```

Numerical interface diffusion or dispersion is a spatial-transport error; it
must not be accompanied by loss of material mass or by spurious changes in the
uniform mechanical and thermal state.

With capillarity absent, a hydrostatic equilibrium satisfies

```{math}
:label: chap13_hydrostatic_balance

\mathbf{u}=\mathbf{0},
\qquad
\nabla p
=
\rho(p,T,\alpha,\widetilde{Y})\mathbf{g}.
```

The pressure gradient therefore changes consistently with the mixture density
through a material interface. A well-balanced discretization evaluates the
pressure gradient and $\rho\mathbf g$ on compatible face/cell geometry so that
their discrete residual vanishes for this state; inserting a cell-centered
density into an otherwise unrelated pressure stencil does not guarantee
hydrostatic balance. With gravity absent and constant pair surface tension
$\sigma_{mn}$, the sharp-interface Young--Laplace equilibrium satisfies

```{math}
:label: chap13_young_laplace_balance

-\nabla p+\mathbf{f}_{\sigma,mn}=\mathbf{0},
\qquad
p_{in}-p_{out}=\sigma_{mn}\kappa_{mn},
```

where material $m$ is inside and material $n$ is outside. With the convention
above, $\mathbf n_{mn}$ points from the outside material toward the inside
material and $\kappa_{mn}>0$ for a convex inner region. The corresponding force
is directed inward. A spherical interface in $d$ spatial dimensions has
$\kappa_{mn}=(d-1)/R$ and
$p_{in}-p_{out}=\sigma_{mn}(d-1)/R$. The discrete pressure and capillary forces
must approach the same balance; otherwise a nominally static interface
develops parasitic motion.

The local conservation laws also impose integral identities. On a fixed
domain $\Omega$ with outward normal $\mathbf{n}$,

```{math}
:label: chap13_material_integral_balance

\frac{d}{dt}\int_{\Omega}q_m\,dV
=
-\int_{\partial\Omega}q_m\mathbf{u}\cdot\mathbf{n}\,dA
+\int_{\Omega}S_m\,dV,
```

and, for $k\in\mathcal{K}_m$,

```{math}
:label: chap13_species_integral_balance

\frac{d}{dt}
\int_{\Omega}q_m\widetilde{Y}^{(m)}_k\,dV
=
-\int_{\partial\Omega}
\left(
q_m\widetilde{Y}^{(m)}_k\mathbf{u}
+\boldsymbol{\mathcal{J}}^{(m)}_k
\right)\cdot\mathbf{n}\,dA
+\int_{\Omega}
\left(\alpha_m\dot\omega^{(m)}_k+S^{tr}_{mk}\right)dV.
```

Consequently, material and species masses are constant in a closed domain
when their corresponding sources vanish, while closed transfer preserves the
summed total mass. In the absence of body forces, capillarity,
and external energy sources, an impermeable, adiabatic, no-work boundary also
gives

```{math}
:label: chap13_total_energy_integral_balance

\frac{d}{dt}\int_{\Omega}(\rho h_t-p)\,dV
=
\frac{d}{dt}\int_{\Omega}\rho E\,dV
=0.
```

If $\mathbf{g}=-\nabla\Phi_g$ for a time-independent potential and constant
pairwise surface tension is retained, the source-free closed-system invariant
includes gravitational potential and interfacial energy for closed material
interfaces with no contact-line work,

```{math}
:label: chap13_mechanical_energy_integral_balance

\frac{d}{dt}
\left[
\int_{\Omega}\left(\rho E+\rho\Phi_g\right)dV
+\sum_{m<n}\sigma_{mn}A_{\Gamma,mn}
\right]
=0,
```

provided their work terms use the same force definitions as momentum.

When $\mathbf g=-\nabla\Phi_g$, $\partial_t\Phi_g=0$, and local continuity
contains $\dot m_{tot}$, the gravitational-work identity becomes

```{math}
:label: chap13_gravity_mass_source_identity

\rho\mathbf g\cdot\mathbf u
=
-\frac{\partial(\rho\Phi_g)}{\partial t}
-\nabla\cdot(\rho\mathbf u\Phi_g)
+\Phi_g\dot m_{tot}.
```

Thus every local mass source exchanges gravitational potential energy with the
flow. This includes donor and receiver sources created by phase-change
redistribution. The corresponding energy source contains
$-\Phi_g\dot m_{tot}$ when fluid plus gravitational potential energy is to be conserved.

## Required Mathematical Invariants

The framework is consistent only if all of the following are enforced:

1. Material admissibility:
   $0\le\alpha_m\le1$ and $\sum_m\alpha_m=1$.
2. Owner-species consistency: the exact material state satisfies
   $\widetilde Y_k^{(m)}\ge0$,
   $\sum_{k\in\mathcal K_m}\widetilde Y_k^{(m)}=1$, and
   $\sum_{k\in\mathcal K_m}z_{mk}=q_m$. Material mass is defined by
   Eq. {eq}`chap13_material_species_mass_sum`, with every species balance
   retained. Global $Y_k$ follows from the material mapping.
3. VOF/species compatibility: Eq. {eq}`chap13_vof_species_manifold`.
4. Appearance consistency: every accepted positive material amount has
   matching species and energy supplied by accepted history, material-bearing
   inflow, or a complete physical source. If all the material's species
   densities are zero, its composition and EOS state are undefined; see
   {ref}`chap13_material_absence_and_appearance`.
5. PTE consistency: every active material EOS is evaluated at the common
   pressure and temperature with its own material composition. At the nonlinear
   fixed point, both residuals in Eq. {eq}`chap13_pte_closure_residuals` vanish;
   Eq. {eq}`chap13_temperature_from_enthalpy` is the fixed-species-mass
   subproblem used to approach that coupled state.
6. Pressure-volume consistency: the pressure-correction equation enforces
   volume closure for the species-derived material masses through
   Eq. {eq}`chap13_fixed_q_volume_error`, applies the fixed-enthalpy response
   in Eq. {eq}`chap13_fixed_q_volume_response`, and converts the
   summed species flux response to each active receiver row's volume weight
   with Eq. {eq}`chap13_fixed_q_receiver_face_weight`. Equation
   {eq}`chap13_fixed_q_pressure_row` then enforces the derived material-volume
   constraint.
7. Transport consistency: Eqs. {eq}`chap13_material_flux`,
   {eq}`chap13_species_flux`,
   {eq}`chap13_species_face_flux_sum`,
   {eq}`chap13_mixture_material_flux_identity`,
   {eq}`chap13_incoming_picard_mass_flux`,
   {eq}`chap13_material_momentum_flux_identity`,
   {eq}`chap13_material_enthalpy_flux_identity`, and
   {eq}`chap13_incoming_direct_z_enthalpy_flux`--{eq}`chap13_material_enthalpy_deferred_correction`.
8. Material-block diffusion conservation and support:
   Eqs. {eq}`chap13_material_diffusion_support` and
   {eq}`chap13_material_species_diffusion_flux_sum` hold on every face.
9. Energy consistency: the same physical species-diffusion flux is reused in
   Eq. {eq}`chap13_cmim_species_enthalpy_face_flux`; the pressure history,
   boundary fluxes, and momentum-force work satisfy
   Eqs. {eq}`chap13_total_enthalpy_bdf1_pressure_time`,
   {eq}`chap13_total_energy_boundary_flux`, and
   {eq}`chap13_static_total_enthalpy_source_relation`.
10. Source consistency: Eqs. {eq}`chap13_material_transfer_source_constraint`,
   {eq}`chap13_phase_change_source_closure`, and
   {eq}`chap13_phase_change_energy_closure`--{eq}`chap13_phase_change_total_energy_source`,
   together with the associated momentum and volume sources.
11. Limiting-state consistency: the stationary and translating contact,
    hydrostatic, Young--Laplace, and integral balances in
    Eqs. {eq}`chap13_stationary_contact_state`--{eq}`chap13_mechanical_energy_integral_balance`.

A numerical scheme must preserve these conservation identities and produce
finite thermodynamic states. FOU and coupled SOU preserve the donor material
sum at each face. MIG reconstructs temporary unit-sum fractions and preserves
that unit sum through its coupled slope balancing. The material-local reconstruction gives zero-sum explicit composition
corrections. Each material convective flux is the sum of its species fluxes,
with the same reconstructed material coefficient in every species row. At
finite Picard iteration, the implicit diffusion corrections need not yet sum
to zero; this requires outer convergence. The converged species state must
be nonnegative, and its EOS-derived volume fractions must satisfy their bounds
and unit-sum constraint.



## Admissible Initial State

An initial condition must define the primary material and thermodynamic state
in each cell:

1. the materials and their species ownership sets,
2. the initial volume fractions $\alpha_m$, with $\sum_m \alpha_m=1$,
3. the initial material compositions $\widetilde{Y}^{(m)}_k$, with
   $\sum_{k\in\mathcal{K}_m}\widetilde{Y}^{(m)}_k=1$ for every active material,
4. the initial pressure and temperature and the assigned EOS branch of every
   active material,
5. the initial velocity field.

The initialization sequence is:

1. build each material composition vector from the specified owner-material
   species list,
2. verify that the common initial pressure and temperature lie in the
   intersection of the active-material EOS domains,
3. evaluate every active material state on its assigned branch and initialize
   $z_{mk}=\alpha_m\rho_m\widetilde Y_k^{(m)}$ for every owned species;
   derive $q_m=\sum_k z_{mk}$, which equals $\alpha_m\rho_m$ because the
   material composition sums to one; set an absent material's species
   densities to zero without querying its composition or EOS,
4. compute the mixture density and material mass fractions,
5. derive the mixture species mass fractions with
   Eq. {eq}`chap13_global_from_internal`,
6. form $e_m$, $e$, $h$, $h_t$, and $E$ from
   Eqs. {eq}`chap13_material_internal_energy`--{eq}`chap13_mixture_total_energy`
   and the initial velocity, and verify that the initialized conserved energy
   satisfies $\rho h_t-p=\rho E$.

This sequence supplies the initial state for the assigned-branch PTE solve.

For an absent material $m_a$, every owned species density is zero, giving
$q_{m_a}=\alpha_{m_a}=0$. Its composition and material EOS state are
undefined. Material appearing through transport receives its species from
the incoming fluxes. Material created by transfer receives its species mass
from $S_{m_a k}^{tr}$ and its energy from the corresponding transfer source.



```{include} material_appearance_face_reference.inc
```

```{include} transport_variants.inc
```


## References

1. S. Thakur, J. Wright, and C. Neal, "An Algebraic Volume of Fluid Method
   Based on NVD and TVD Schemes Coupled to an All-Speed CFD Solver," AIAA
   SciTech Forum (2022), doi:10.2514/6.2022-1203.
2. J. U. Brackbill, D. B. Kothe, and C. Zemach, "A continuum method for
   modeling surface tension," J. Comput. Phys. 100, 335-354 (1992).
3. M. Francois et al., "A Balanced-Force Algorithm for Continuous and Sharp
   Interfacial Surface Tension Models Within a Volume Tracking Framework," J.
   Comput. Phys. 213, 141-173 (2006).
4. F. Denner and B. van Wachem, "Fully-Coupled Balanced-Force VOF Framework
   for Arbitrary Meshes with Least-Squares Curvature Evaluation from Volume
   Fractions," Numer. Heat Transfer B 65, 218-255 (2014).
5. A. D. Demou, N. Scapin, M. Pelanti, and L. Brandt, "A pressure-based
   diffuse interface method for low-Mach multiphase flows with mass transfer,"
   J. Comput. Phys. 448, 110730 (2022), doi:10.1016/j.jcp.2021.110730.
6. B. Clayton, J. McConnell, and Clell J. Solomon, "Analysis, thermodynamics, and a
   numeric solver for a pressure-temperature equilibrium closure of the
   four-equation model," arXiv:2606.27726 (2026).
7. C. Kunkelmann, *Numerical Modeling and Investigation of Boiling
   Phenomena*, doctoral dissertation, Technische Universität Darmstadt (2011),
   doi:10.26083/tuprints-00002731.
8. S. Hardt and F. Wondra, "Evaporation model for interfacial flows based on a
   continuum-field representation of the source terms," J. Comput. Phys. 227,
   5871--5895 (2008), doi:10.1016/j.jcp.2008.02.020.
9. P. J. Zwart, A. G. Gerber, and T. Belamri, "A Two-Phase Flow Model for
   Predicting Cavitation Dynamics," Fifth International Conference on
   Multiphase Flow, Paper 152 (2004).
10. R. W. Schrage, *A Theoretical Study of Interphase Mass Transfer*, Columbia
    University Press (1953).