15. 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:

One common pressure, temperature, and velocity in each cell. Material closures are evaluated at this common $(p,T)$ and with their own material compositions.
Material interfaces are tracked by VOF fields. Passive species diffusion does not transfer species between different materials.
Species belong to exactly one material set and mix only with species owned by that same material.
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.
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.

15.1. 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:

$m = 1,...,N_m$ (index of materials)
$k = 1,...,N_s$ (index of species), with $N_{s,m}=|\mathcal K_m|$ species owned by material $m$
$\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:

(15.1.1)#\[\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:

(15.1.2)#\[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

(15.1.3)#\[\begin{split}\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}\end{split}\]

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:

(15.1.4)#\[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

(15.1.5)#\[\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

(15.1.6)#\[\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

(15.1.7)#\[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

(15.1.8)#\[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

(15.1.9)#\[Y_k=w_{m(k)}\widetilde{Y}^{(m(k))}_k\]

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

(15.1.10)#\[\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:

(15.1.11)#\[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\]

Meaning of simplex in this chapter

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

(15.1.12)#\[\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.

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

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

(15.2.1.1)#\[\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:

(15.2.1.2)#\[\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$.

15.2.1.1. Outer iteration#

The general pressure-based iteration is defined in Chapter 3. 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 Chapter 4:

(15.2.1.3)#\[{}^{(\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 Experimental material-to-energy predictor coupling changes this dependency; it is separate from the standard sequence described here.

The information flow through one CMIM outer iteration is

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 Segregated Nonlinear Algorithm.

15.2.2. 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:

(15.2.2.1)#\[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. (15.5.1.6).

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. (15.2.1.2). 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 Material-Appearance Face Reference.

15.2.3. 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. (15.3.2)
Material enthalpy	$h_m$	Eq. (15.3.3)
Material internal energy	$e_m$	Eq. (15.3.8)
Mixture density	$\rho$	Eq. (15.3.9)
Mixture static thermodynamic enthalpy	$h$	Eq. (15.3.10)
Mixture internal energy	$e$	Eq. (15.3.11)
Mixture total energy	$E$	Eq. (15.3.12)
Mixture species mass fraction	$Y_k$	Eq. (15.1.7)
Material mass fraction	$w_m$	Eq. (15.1.8)
Fixed-enthalpy volume compressibility	$\chi$	Eq. (15.7.1.3)

15.3. 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:

(15.3.1)#\[\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:

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

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

It also supplies the fixed-composition derivatives

(15.3.4)#\[\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

(15.3.5)#\[\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

(15.3.6)#\[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:

(15.3.7)#\[\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

(15.3.8)#\[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:

(15.3.9)#\[\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:

(15.3.10)#\[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:

(15.3.11)#\[\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

(15.3.12)#\[E=e+K, \qquad \rho h_t-p=\rho E.\]

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

(15.3.13)#\[h_t=h+\frac{1}{2}\lVert\vec{u}\rVert^2\]

15.3.1. 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. (15.2.1.2). 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

(15.3.1.1)#\[\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 (15.3.10) also requires $\sum_{m\in\mathcal A}q_mh_m=\rho h^{\mathrm{tar}}$. These two requirements define the closure residuals

(15.3.1.2)#\[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

(15.3.1.3)#\[\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

(15.3.1.4)#\[\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

(15.3.1.5)#\[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. (15.7.1.3). 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 Segregated Nonlinear Algorithm.

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

(15.4.1)#\[\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

(15.4.2)#\[\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

(15.4.3)#\[\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. (4.2.1.16). The two-point conductive energy rate, positive from $L$ to $R$, is

(15.4.4)#\[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. (15.9.6.4).

15.5. Governing Equations#

15.5.1. 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:

(15.5.1.1)#\[\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. (15.5.1.1) is

(15.5.1.2)#\[\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$.

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

(15.5.1.3)#\[\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

(15.5.1.4)#\[\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$,

(15.5.1.5)#\[\sum_{k\in\mathcal{K}_m}\dot{\omega}^{(m)}_k=0.\]

The transfer sources satisfy

(15.5.1.6)#\[\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}.\]

15.5.2. Material and mixture mass balances#

Summing Eq. (15.5.1.1) over the species of material $m$ gives

(15.5.2.1)#\[\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

(15.5.2.2)#\[\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:

(15.5.2.3)#\[\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.

15.5.3. Mixture momentum#

(15.5.3.1)#\[\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.

15.5.4. Mixture static enthalpy#

Define the material-resolved species-enthalpy diffusion flux

(15.5.4.1)#\[\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

(15.5.4.2)#\[\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. (15.5.3.1), where $\mathbf{S}_u$ is the additional conservative momentum source and gravity and capillarity are shown separately, the coupled source terms satisfy

(15.5.4.3)#\[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. (15.3.1.5), 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. (15.1.7) does not replace the independently evaluated material thermodynamic states.

15.5.5. Mixture total enthalpy#

The mixture total enthalpy is defined by Eq. (15.3.13). Its conservative continuum balance is

(15.5.5.1)#\[\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:

(15.5.5.2)#\[\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. (15.5.5.1) is

(15.5.5.3)#\[\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}$.

15.5.5.1. Discrete enthalpy flux due to species diffusion#

At an internal face, the energy equation reuses the species flux from Eq. (15.6.6.14):

(15.5.5.4)#\[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

(15.5.5.5)#\[\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:

(15.5.5.6)#\[\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.

15.5.6. Constitutive and term definitions#

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

(15.5.6.1)#\[\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. (15.5.4.2) is:

(15.5.6.2)#\[\Phi=\boldsymbol{\tau}:\nabla\mathbf{u}\]

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

(15.5.7.1)#\[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

(15.5.7.2)#\[\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. (15.5.7.2) 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

(15.5.7.3)#\[\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

(15.5.7.4)#\[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

(15.5.7.5)#\[\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

(15.5.7.6)#\[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. (15.5.7.6) supplies the preliminary curvature

(15.5.7.7)#\[\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

(15.5.7.8)#\[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

(15.5.7.9)#\[\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

(15.5.7.10)#\[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

(15.5.7.11)#\[\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

(15.5.7.12)#\[\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,

(15.5.7.13)#\[\mathbf n_{mn}\delta_{\Gamma,mn}=\nabla C_{mn}.\]

The surface-tension force density is therefore

(15.5.7.14)#\[\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

(15.5.7.15)#\[\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.

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

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

(15.5.7.16)#\[\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,

(15.5.7.17)#\[\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.

15.5.8. Source terms and sign conventions#

Term definitions and sign conventions used in Eqs. (15.5.2.3), (15.5.2.2), (15.5.1.1), (15.5.3.1), (15.5.4.2), and (15.5.5.1):

$\dot{m}_m$ is the net material mass source per mixture volume for material $m$; positive values add material mass.
$\dot{\omega}^{(m)}_k$ is an intrinsic species-production rate per material volume and is multiplied by $\alpha_m$.
$S^{tr}_{mk}$ is a conservative inter-material or external species source per mixture volume and is not multiplied by $\alpha_m$ again.
$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.
$\mathbf{f}_{\sigma}$ is volumetric surface-tension force density.
$\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.

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

15.6.1. Face fluxes and cell balances#

Following the solver framework in Chapter 3 and the finite-volume derivation in Chapter 4, 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:

(15.6.1.1)#\[F_{z,ms,f}=\dot{\Omega}_f z_{ms,f}\]

Here $\dot{\Omega}_f$ is the face volume flow rate defined in Eq. (5.1), with units of volume per time. Consequently $F_{z,ms,f}$ has units of mass per time. The generic scalar equation in Eq. (4.2.3) 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:

(15.6.1.2)#\[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 Loci Face Orientation and Cell Incidence converts these stored fluxes to the outward orientation of cell $P$:

(15.6.1.3)#\[\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. (15.5.1.1), over the fixed cell and applying BDF1 gives

(15.6.1.4)#\[\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:

(15.6.1.5)#\[\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

(15.6.1.6)#\[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.

15.6.2. Volume-fraction face reconstruction#

15.6.2.1. 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 (15.6.1.3) supplies its sign in the balance for cell $P$. First-order upwinding (FOU) uses

(15.6.2.1)#\[\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

(15.6.2.2)#\[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:

(15.6.2.3)#\[\begin{split}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}\end{split}\]

The raw SOU candidate is

(15.6.2.4)#\[\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.

15.6.2.2. 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:

(15.6.2.5)#\[\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 Consistent material fractions for MIG 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

(15.6.2.6)#\[\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

(15.6.2.7)#\[\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

(15.6.2.8)#\[\begin{split}\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}\end{split}\]

and form the face gradient by opposite-volume weighting,

(15.6.2.9)#\[\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,

(15.6.2.10)#\[\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$.

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

(15.6.2.11)#\[a_{m,D}^{SOU}=\alpha_{m,D} \qquad a_{m,D}^{MIG}=\beta_{m,D}\]

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

(15.6.2.12)#\[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:

(15.6.2.13)#\[\begin{split}\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}\end{split}\]

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

(15.6.2.14)#\[\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,

(15.6.2.15)#\[\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.

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

(15.6.2.16)#\[\begin{split}\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}\end{split}\]

Composition uses the same support condition. Its SOU increments are projected to zero sum within each material block, as defined in Eq. (15.6.4.5). 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$.

15.6.3. 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:

(15.6.3.1)#\[\begin{split}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}\end{split}\]

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.

15.6.4. Material composition at a face#

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

(15.6.4.1)#\[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:

(15.6.4.2)#\[\begin{split}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}\end{split}\]

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

(15.6.4.3)#\[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:

(15.6.4.4)#\[\begin{split}\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}\end{split}\]

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

(15.6.4.5)#\[\begin{split}\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.\end{split}\]

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

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

(15.6.5.1)#\[\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. (15.6.4.5). Multiplying this composition by the reconstructed material density gives the reconstructed species density and its advective flux:

(15.6.5.2)#\[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}\]

15.6.5.1. From the reconstructed flux to the standard predictor#

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

(15.6.5.3)#\[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

(15.6.5.4)#\[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. (4.2.2.12), 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:

(15.6.5.5)#\[\begin{split}\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}\end{split}\]

In the species balance, Eq. (15.6.1.4), 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.

15.6.5.2. 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}$:

(15.6.5.6)#\[\begin{split}\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}\end{split}\]

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

(15.6.5.7)#\[\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:

(15.6.5.8)#\[\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

(15.6.5.9)#\[\begin{split}\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}\end{split}\]

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

(15.6.5.10)#\[\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

(15.6.5.11)#\[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:

(15.6.5.12)#\[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.

15.6.5.3. 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 Implicit material advection. The diffusion and cell-balance construction below apply to either convection split.

Conservative rescaling, developed in Conservative flux 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. (15.15.2.12). The flux and pressure relations below use these selected quantities.

15.6.6. Finite-volume species diffusion#

The generic finite-volume diffusion split is developed in Decomposition of Generic Transport Equation to Algebraic Cell Equation. 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

(15.6.6.1)#\[\Gamma_{mk}=\alpha_m\rho_mD_{mk}.\]

A material face link is supported only when

(15.6.6.2)#\[\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

(15.6.6.3)#\[\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

(15.6.6.4)#\[\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

(15.6.6.5)#\[\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

(15.6.6.6)#\[\begin{split}\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}\end{split}\]

The positive compact conductance is therefore

(15.6.6.7)#\[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

(15.6.6.8)#\[\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

(15.6.6.9)#\[\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. (15.6.6.5) into the integrated Fick flux:

(15.6.6.10)#\[\begin{split}\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}\end{split}\]

The complete raw integrated face flux is

(15.6.6.11)#\[\begin{split}\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}\end{split}\]

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

(15.6.6.12)#\[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

(15.6.6.13)#\[\begin{split}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}\end{split}\]

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

(15.6.6.14)#\[\begin{split}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}\end{split}\]

The sign-group correction is applied once to the complete raw flux in Eq. (15.6.6.11), 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

(15.6.6.15)#\[\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

(15.6.6.16)#\[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

(15.6.6.17)#\[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

(15.6.6.18)#\[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. (15.6.6.17) is positive from $L$ to $R$. Its cell-outward form is

(15.6.6.19)#\[\begin{split}F_{ms,P,f}^{lin,*,k} = \begin{cases} +F_{ms,f}^{lin,*,k}, & P=L,\\ -F_{ms,f}^{lin,*,k}, & P=R. \end{cases}\end{split}\]

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. (15.6.6.15); 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.

15.6.6.1. 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. (15.6.6.11) 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.

15.6.7. Species predictor row#

Use the standard convective flux in Eq. (15.6.5.4), or the optional implicit material flux in Eq. (15.15.1.9), together with the linearized diffusion flux in Eq. (15.6.6.17). The predictor form of the BDF1 species balance, Eq. (15.6.1.4), is

(15.6.7.1)#\[{}^{(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:

(15.6.7.2)#\[{}^{(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.

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

(15.6.7.3)#\[\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. (15.6.1.5).

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:

(15.6.7.4)#\[\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.

15.7. 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. (15.3.1.5) 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

(15.7.1)#\[\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.

15.7.1. Predicted material-volume closure error#

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

(15.7.1.1)#\[\begin{split}E_{V,P}^{\mathrm{pred},k} = \sum_{\substack{m\\q_{m,P}^{*,k}>0}} \alpha_{m,P}^{\mathrm{pred},k} -1.\end{split}\]

$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. (15.7.1.1) 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

(15.7.1.2)#\[\begin{split}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}.\end{split}\]

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

(15.7.1.3)#\[\begin{split}\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].\end{split}\]

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

(15.7.1.4)#\[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.

15.7.2. Face-flux response#

With implicit material advection, the material convective flux is the sum of the species fluxes in Eq. (15.15.1.10). Its explicit composition correction has zero material sum. Holding the donor, reconstructed coefficient, and predicted species densities fixed gives

(15.7.2.1)#\[\begin{split}\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}\end{split}\]

This response applies to the implicit material linearization. For standard deferred material advection, use Eq. (15.15.1.11). 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

(15.7.2.2)#\[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:

(15.7.2.3)#\[\begin{split}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}.\end{split}\]

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

(15.7.2.4)#\[\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

(15.7.2.5)#\[\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.

15.7.3. 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:

(15.7.3.1)#\[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

CMIM Pressure-Correction Row

(15.7.3.2)#\[\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.

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

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

(15.7.4.1)#\[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.

(15.7.4.2)#\[(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:

(15.7.4.3)#\[\begin{split}\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}\end{split}\]

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

(15.7.4.4)#\[\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.

15.8. 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. (15.6.1.2). 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. (15.6.4.5) are applied separately to each material.

Six finite-volume cells showing a notional helium-nitrogen boundary, an internal face, and upward material transport. — Fig. 15.8.1 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. (15.6.2.16) and (15.6.4.5). The curved helium–nitrogen boundary visualizes the cell-average occupancies; it does not represent an explicitly reconstructed geometric interface.#

15.8.1. Material and species fluxes#

The material mass flux associated with the reconstructed state is

(15.8.1.1)#\[\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. (15.15.2.5). The species fluxes are

(15.8.1.2)#\[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. (15.15.2.15). Consequently,

(15.8.1.3)#\[\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. (15.15.1.10) 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.

15.8.2. Mixture mass and momentum fluxes#

The mixture mass flux is the sum over materials:

(15.8.2.1)#\[\dot m_f=\sum_m\dot m_{m,f}.\]

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

(15.8.2.2)#\[\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

(15.8.2.3)#\[\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.

15.8.3. 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. (15.15.1.9) at the incoming state gives

(15.8.3.1)#\[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

(15.8.3.2)#\[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:

(15.8.3.3)#\[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

(15.8.3.4)#\[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

(15.8.3.5)#\[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

(15.8.3.6)#\[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

(15.8.3.7)#\[\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. (15.5.5.4).

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

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

(15.9.1.1)#\[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,

(15.9.1.2)#\[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,

(15.9.1.3)#\[\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

(15.9.1.4)#\[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

(15.9.1.5)#\[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

(15.9.1.6)#\[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

(15.9.1.7)#\[(\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

(15.9.1.8)#\[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.

15.9.2. 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. (15.9.1.8). For either side $s\in\{d,r\}$, write the heat entering the interface as

(15.9.2.1)#\[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:

(15.9.2.2)#\[\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

(15.9.2.3)#\[\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:

(15.9.2.4)#\[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$:

(15.9.2.5)#\[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. (15.9.1.5), 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. (15.9.2.4) 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.

15.9.3. 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$:

(15.9.3.1)#\[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. (15.9.3.1) 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:

(15.9.3.2)#\[\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

(15.9.3.3)#\[\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.

15.9.4. 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:

(15.9.4.1)#\[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.

15.9.5. Hardt–Wondra conservative redistribution#

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

(15.9.5.1)#\[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:

(15.9.5.2)#\[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:

(15.9.5.3)#\[\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

(15.9.5.4)#\[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:

(15.9.5.5)#\[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:

(15.9.5.6)#\[\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.

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

(15.9.6.1)#\[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

(15.9.6.2)#\[\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:

(15.9.6.3)#\[\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

(15.9.6.4)#\[\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. (15.9.3.2). 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

(15.9.6.5)#\[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. (15.5.4.3) 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:

(15.9.6.6)#\[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.

15.9.7. 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,

(15.9.7.1)#\[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,

(15.9.7.2)#\[ \begin{align}\begin{aligned}\dot m_c^{vol} = \begin{cases}\\\begin{split}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]\end{split}\\\begin{split}-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]\end{split}\\0, & p=p_{sat}(T).\\\end{cases}\end{aligned}\end{align} \]

$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. (15.9.1.6). 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.

15.10. Segregated Nonlinear Algorithm#

Each outer iteration follows the pressure-based procedure in 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.

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$.
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. (15.9.4.1) fixed during the nonlinear solve.
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.
Solve the species, momentum, and total-enthalpy predictors. The momentum and energy equations use the incoming mass flux in Eq. (15.8.2.2); the enthalpy correction uses Eq. (15.8.3.2). No predictor requires another predictor’s newly solved value.
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. (15.3.1.5) for $T^{\mathrm{pred},k}$. Evaluate the predicted EOS densities, material volumes, and fixed-enthalpy volume response.
Use the predicted volume error and the difference between the momentum and species predictor volume fluxes to assemble Eq. (15.7.3.2). Solve for $p^{\prime,k}$.
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}$.
Check conservation-equation residuals, the species update in Eq. (15.6.7.4), material-volume closure, enthalpy closure, and continuity using Eq. (15.7.4.4). Account for pressure-induced material appearance using Eq. (15.7.2.4). 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.

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

(15.11.1)#\[\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

(15.11.2)#\[\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

(15.11.3)#\[\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

(15.11.4)#\[-\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}$,

(15.11.5)#\[\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$,

(15.11.6)#\[\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

(15.11.7)#\[\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,

(15.11.8)#\[\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

(15.11.9)#\[\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.

15.12. Required Mathematical Invariants#

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

Material admissibility: $0\le\alpha_m\le1$ and $\sum_m\alpha_m=1$.
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. (15.1.4), with every species balance retained. Global $Y_k$ follows from the material mapping.
VOF/species compatibility: Eq. (15.1.10).
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 Material absence and appearance.
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. (15.3.1.2) vanish; Eq. (15.3.1.5) is the fixed-species-mass subproblem used to approach that coupled state.
Pressure-volume consistency: the pressure-correction equation enforces volume closure for the species-derived material masses through Eq. (15.7.1.1), applies the fixed-enthalpy response in Eq. (15.7.1.3), and converts the summed species flux response to each active receiver row’s volume weight with Eq. (15.7.2.3). Equation (15.7.3.2) then enforces the derived material-volume constraint.
Transport consistency: Eqs. (15.8.1.1), (15.8.1.2), (15.8.1.3), (15.8.2.1), (15.8.2.2), (15.8.2.3), (15.8.3.4), and (15.8.3.2)–(15.8.3.7).
Material-block diffusion conservation and support: Eqs. (15.6.6.2) and (15.6.6.15) hold on every face.
Energy consistency: the same physical species-diffusion flux is reused in Eq. (15.5.5.4); the pressure history, boundary fluxes, and momentum-force work satisfy Eqs. (15.5.5.2), (15.5.5.3), and (15.5.4.3).
Source consistency: Eqs. (15.5.1.6), (15.9.5.6), and (15.9.6.4)–(15.9.6.5), together with the associated momentum and volume sources.
Limiting-state consistency: the stationary and translating contact, hydrostatic, Young–Laplace, and integral balances in Eqs. (15.11.1)–(15.11.8).

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.

15.13. Admissible Initial State#

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

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

The initialization sequence is:

build each material composition vector from the specified owner-material species list,
verify that the common initial pressure and temperature lie in the intersection of the active-material EOS domains,
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,
compute the mixture density and material mass fractions,
derive the mixture species mass fractions with Eq. (15.1.7),
form $e_m$, $e$, $h$, $h_t$, and $E$ from Eqs. (15.3.8)–(15.3.12) 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.

15.14. Material-Appearance Face Reference#

This section applies Material absence and appearance to one material $m$ at an internal face. An empty cell means empty of material $m$; other materials still occupy the cell. Local material presence does not by itself establish that a wider reconstruction stencil supports that material.

15.14.1. One oriented internal face#

Let $P=L_f$ and $N_f=R_f$, with $\mathbf S_f$ pointing from $P$ to $N_f$. For all four sketches, take $\dot{\Omega}_f^k<0$. The donor is then $D_f^k=N_f$, and a negative material mass flux is inflow to $P$:

(15.14.1.1)#\[\dot{\Omega}_f^k=\mathbf u_f^k\mathbin{\cdot}\mathbf S_f, \qquad \dot m_{m,f}^k=\widehat q_{m,f}^k\dot{\Omega}_f^k.\]

Here $\widehat q_{m,f}^k$ is the selected material face density. In standard mode it equals the reconstructed $\alpha_{m,f}^k\rho_{m,f}^k$ for a present donor and is zero for an absent donor.

15.14.1.1. Material absent from both cells#

$Adjacent cells P and N f share face f. Cell-center dots mark P and N f, and material m has zero volume fraction in both cells.$

Fig. 15.14.1.1 Material $m$ is absent at both cell centers.#

The donor is empty, so $\alpha_{m,f}^k=\widehat q_{m,f}^k=0$. The incoming state gives no advective material, species, or enthalpy flux for $m$ at this face. History, sources, or another face may still supply species during the implicit solve.

15.14.1.2. Material present only in the receiver#

$Adjacent cells P and N f share face f. Cell-center dots mark P and N f. Material m has positive volume fraction in P and zero volume fraction in N f.$

Fig. 15.14.1.2 Material $m$ is present in receiver $P$ and absent from donor $N_f$.#

The empty donor again gives $\alpha_{m,f}^k=\widehat q_{m,f}^k=0$. This face supplies no incoming-state inflow of $m$. Material already in $P$, its accepted history, and its local sources remain in the cell balance. Reversing the volume flux would make $P$ the donor.

15.14.1.3. Material present in both cells#

$Adjacent cells P and N f share face f. Cell-center dots mark P and N f, and material m has positive volume fraction in both cells.$

Fig. 15.14.1.3 Material $m$ is present at both cell centers.#

With FOU reconstruction, $\widehat q_{m,f}^k=q_{m,N_f}^k$. A higher-order VOF scheme may give a different material face density. SOU density, composition, and species enthalpy require material support throughout the donor’s gradient stencil, not just in these two cells. For $\widehat q_{m,f}^k>0$, the negative volume flux carries material into $P$.

15.14.1.4. Material present only in the donor#

$Adjacent cells P and N f share face f. Cell-center dots mark P and N f. Material m has zero volume fraction in P and positive volume fraction in N f.$

Fig. 15.14.1.4 Material $m$ is absent from receiver $P$ and present in donor $N_f$.#

This is the appearance case. The VOF scheme supplies the material face fraction. Density, composition, and species enthalpies use donor values because the conditional-state stencil crosses into a cell where material $m$ is absent. The species equations carry these donor species into $P$ through their common material coefficient.

No receiver composition or intrinsic density is needed by this advective flux. Once the species predictor has a positive material sum, its composition and the common pressure and temperature define the new receiver EOS state.

Reversing $\dot{\Omega}_f^k$ interchanges donor and receiver. At zero volume flux, the incoming-state advective flux is zero.

15.14.2. Lagged face states and the implicit species solve#

An absent lagged donor does not force the implicit predictor flux to remain zero throughout the global solve. For $q_{m,D_f^k}^k=0$, Eq. (15.15.1.9) gives

(15.14.2.1)#\[F_{z,ms,f}^{*,k}=\dot{\Omega}_f^kz_{ms,D_f^k}^{*,k}, \qquad F_{q,m,f}^{*,k} =\sum_{s\in\mathcal K_m}F_{z,ms,f}^{*,k} =\dot{\Omega}_f^kq_{m,D_f^k}^{*,k}.\]

Another incoming face can make the donor species predictor positive. The same implicit coupling can then transport species through this face and through several cells that were empty in the incoming iterate. Material mass follows from the species sum at every stage.

The reconstructed coefficients and enthalpy fluxes remain lagged during the predictor solves and are rebuilt at the next outer iteration. At convergence, incoming and predicted states agree. A donor that remains exactly empty then has zero material, species, and enthalpy flux.

15.14.3. The appearing material in the pressure row#

After the species solve, derive $q_{m,P}^{*,k}=\sum_s z_{ms,P}^{*,k}$. If this sum is zero, the material does not contribute to the cell’s predicted volume or pressure response. If it is positive, its composition is $z_{ms,P}^{*,k}/q_{m,P}^{*,k}$, and the pre-pressure closure supplies $\rho_{m,P}^{\mathrm{pred},k}$.

The pressure response follows the selected species linearization. With implicit material advection and a present lagged donor,

(15.14.3.1)#\[\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f} = \frac{\widehat q_{m,f}^k}{q_{m,N_f}^k}q_{m,N_f}^{*,k}.\]

For an absent lagged donor it is $q_{m,N_f}^{*,k}$, as in Eq. (15.7.2.1). Standard deferred material advection uses $q_{m,N_f}^{*,k}$ for either donor state, following Eq. (15.15.1.11). The row’s volume-flux weight is therefore

(15.14.3.2)#\[\begin{split}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}.\end{split}\]

The receiver density converts the material mass-flux response to a volume-flux response. The cell source is the predicted volume error

(15.14.3.3)#\[\begin{split}E_{V,P}^{\mathrm{pred},k} = \sum_{\substack{m\\q_{m,P}^{*,k}>0}} \frac{q_{m,P}^{*,k}}{\rho_{m,P}^{\mathrm{pred},k}}-1.\end{split}\]

Both use the species-derived material masses and their predicted EOS states. The species densities, material reconstruction coefficient, and explicit composition correction stay fixed during pressure correction. The cell thermodynamic response is the fixed-enthalpy derivative in Eq. (15.7.1.3).

If the corrected volume flux activates or reverses a face into an absent receiver material, the appearance estimate in Eq. (15.7.2.4) measures the additional volume change. A change above the update tolerance requires another outer iteration, so species transport and thermodynamic closure account for the new inflow.

15.14.4. Composition supplied by inflow#

Consider a converged, source-free, diffusion-free FOU state with no accepted history of material $m$ in cell $P$. Its incoming species fluxes give

(15.14.4.1)#\[\widetilde{\mathbf Y}_{P}^{(m)} = \frac{ \sum_f F_{q,m,Pf}^{\mathrm{in}}\widetilde{\mathbf Y}_{D_f}^{(m)} }{ \sum_f F_{q,m,Pf}^{\mathrm{in}} }, \qquad F_{q,m,Pf}^{\mathrm{in}} =\max\!\left(-\sigma_{Pf}F_{q,m,f},0\right),\]

provided the denominator is positive. One inflow supplies its donor composition; several inflows supply their material-mass-flux-weighted composition. This follows by summing the species rows and dividing each species balance by their material sum. Accepted history adds its own mass weight. A source that creates material supplies the species and energy contributions specified by its physical model.

During the iteration, the composition is always derived from the local species predictor and its sum. It is not assigned from a selected neighbor.

15.14.5. Operator behavior near material absence#

Operation	Treatment
Species advection	Use the common material coefficient in every owned species row. Retain implicit FOU coupling for an absent lagged donor.
Material face reconstruction	Reconstruct the derived volume fractions. FOU, guarded SOU, and MIG leave an exactly absent donor component at zero.
Density and species-enthalpy SOU	Extrapolate only with material present throughout the donor’s gradient stencil; otherwise use the donor value.
Composition SOU	Apply the same support condition and preserve zero-sum increments separately within each material.
Species diffusion	Require material in both adjacent cells. Retain the compact term but omit the nonorthogonal correction if the wider gradient stencil is unsupported.
Pressure-volume closure	Use positive species sums and their predicted EOS states. An absent material contributes zero volume.
Pressure face response	Convert the summed species mass-flux response to cell volume using the active predicted cell density.

In an admissible state, exact absence means that every species density of the material is zero. A positive material sum instead defines a composition and requires an EOS state, however small that sum is.

15.15. Transport Formulations and Experimental Extensions#

The species equations advance $z_{ms}$, their sums give $q_m$, and the energy equation advances mixture total enthalpy $h_t$. In the segregated iteration of Outer iteration, each predictor uses coefficients and face reconstructions from the same incoming state. The predictors become mutually consistent as the outer iterations converge.

Four changes to this process are considered here. They address different questions:

Approach	Question addressed	What changes
Implicit material advection	How should a reconstructed material flux enter the species equation?	Part of the explicit correction becomes a coefficient of the unknown species density
Conservative flux rescaling	How much of a face correction can be applied without exhausting a cell’s available mass?	Material and composition corrections are reduced by shared face factors
Experimental material-to-energy coupling	How should the energy predictor account for a change in material proportions?	The completed material predictor supplies a correction to the energy update
Consistent MIG reconstruction	How should an incompletely converged material-volume sum affect the reconstructed interface?	Temporary fractions and gradients satisfy the material-volume identities before reconstruction

The first two are available transport formulations. Consistent reconstruction is part of the current MIG scheme. The material-to-energy coupling remains experimental and changes the independence of the predictors; MIG reconstruction retains their incoming Picard state. The notation retains superscript $n$ for the accepted time level, $k$ for the incoming outer iterate, and $*,k$ for the predictor result. The species index is $s\in\mathcal K_m$.

15.15.1. Implicit material advection#

Consider the convective flux of one species of one material through an internal face. Write $F_f$ for its mass flow rate and $z$ for its mass per mixture volume, suppressing the material and species indices in this derivation. The donor $D=D_f^k$ is selected by the incoming volume flow $\dot{\Omega}_f^k$. Superscript $k$ marks known values from the incoming iteration; $*,k$ marks the unknown predictor. The standard deferred flux from Eq. (15.6.5.4) is

(15.15.1.1)#\[F_f^{\mathrm{std},*,k} =\dot{\Omega}_f^k \left[z_D^{*,k}+\left(z_f^k-z_D^k\right)\right]\]

Here $z_f^k$ is the reconstructed face value. The donor term is implicit, while the correction $z_f^k-z_D^k$ stays fixed during the linear solve. If that correction removes more species than a cell’s history and inflow can supply, the solved cell value can become negative. Implicit material advection changes the treatment of the material part of this correction.

15.15.1.1. Separate material transport from composition transport#

For a material present in the donor, the species density is $z=qY$. Here $q=\alpha\rho$ is material mass per mixture volume, with $\rho$ the intrinsic material density. The species mass fraction within that material is $Y$, abbreviating $\widetilde Y_s^{(m)}$. At the face, $q_f^k$ denotes the reconstructed material density $q_{m,f}^{\mathrm{rec},k}$, and $Y_f^k$ denotes the reconstructed composition $\widehat Y_{ms,f}^k$. Thus $z_f^k=q_f^kY_f^k$ and $z_D^k=q_D^kY_D^k$. Adding and subtracting $q_f^kY_D^k$ splits the face correction into two parts:

(15.15.1.2)#\[\begin{split}\begin{aligned} z_f^k-z_D^k &=q_f^kY_f^k-q_D^kY_D^k\\ &=\left(q_f^k-q_D^k\right)Y_D^k +q_f^k\left(Y_f^k-Y_D^k\right) \end{aligned}\end{split}\]

The first part changes how much material crosses the face at the donor composition. The second changes the species proportions of that material flow. Substituting both parts back into the full flux gives

(15.15.1.3)#\[F_f^{\mathrm{std},*,k} =\dot{\Omega}_f^k \left[ z_D^{*,k} +\underbrace{\left(q_f^k-q_D^k\right)Y_D^k}_{\text{material correction}} +\underbrace{q_f^k\left(Y_f^k-Y_D^k\right)}_{\text{composition correction}} \right]\]

This is still the standard deferred flux. Both corrections are known during the solve. The same decomposition is given with full species and material indices in Eq. (15.6.5.6).

15.15.1.2. Make the material correction implicit#

Assume the donor contains the material, so $q_D^k>0$. In the material correction, use $Y_D^k=z_D^k/q_D^k$:

(15.15.1.4)#\[\left(q_f^k-q_D^k\right)Y_D^k =\frac{q_f^k-q_D^k}{q_D^k}z_D^k\]

This equality only rewrites the known material correction. The alternative linearization now replaces its final $z_D^k$ by the unknown $z_D^{*,k}$. The ratio stays fixed, and the composition correction is unchanged. This replacement is a numerical choice, not an algebraic identity. The full modified flux is

(15.15.1.5)#\[F_f^{\mathrm{imp},*,k} =\dot{\Omega}_f^k \left[ z_D^{*,k} +\frac{q_f^k-q_D^k}{q_D^k}z_D^{*,k} +q_f^k\left(Y_f^k-Y_D^k\right) \right]\]

The first two terms now contain the same unknown. Their coefficients combine as

(15.15.1.6)#\[1+\frac{q_f^k-q_D^k}{q_D^k} =\frac{q_f^k}{q_D^k}\]

The resulting implicit material-advection flux is therefore

(15.15.1.7)#\[F_f^{\mathrm{imp},*,k} =\dot{\Omega}_f^k\frac{q_f^k}{q_D^k}z_D^{*,k} +\dot{\Omega}_f^k q_f^k\left(Y_f^k-Y_D^k\right)\]

The material contribution now responds to the unknown species density during the solve. The composition correction remains explicit. The ratio $q_f^k/q_D^k$ uses incoming values and is held fixed, as are the other Picard coefficients. It is common to all species of the material.

15.15.1.3. Checks on the resulting flux#

With FOU material reconstruction, $q_f^k=q_D^k$. The coefficient ratio is one, the material correction is zero, and the standard and implicit formulations coincide at every iteration.

At nonlinear convergence, $z_D^{*,k}=z_D^k=q_D^kY_D^k$. Substituting into the final flux gives

(15.15.1.8)#\[F_f^{\mathrm{imp}} =\dot{\Omega}_f q_fY_D +\dot{\Omega}_f q_f\left(Y_f-Y_D\right) =\dot{\Omega}_f q_fY_f =F_f^{\mathrm{std}}\]

The iteration superscripts have been dropped in this converged expression. Both formulations recover the same reconstructed species flux; they differ in how the linear solves approach it.

For nonnegative $q_f^k$, the modified coefficient retains the upwind sign. With FOU composition, $Y_f^k=Y_D^k$, so the remaining explicit correction also vanishes. SOU composition can still supply a correction, and diffusion and sources retain their own admissibility requirements. The change removes the deferred material contribution as a source of negative predictors, but is not a general positivity guarantee for the species equation.

15.15.1.4. Selected face density and absent donors#

Restore the material and species indices for the following relations. The explicit composition term in the derived flux is $\delta F_{z,ms,f}^k$ from Eq. (15.6.5.9): $\dot{\Omega}_f^k q_{m,f}^{\mathrm{rec},k}\delta\widetilde Y_{ms,f}^k$, where $\delta\widetilde Y_{ms,f}^k$ is the face-minus-donor composition increment. The material density denoted by $q_f^k$ in the derivation is $q_{m,f}^{\mathrm{rec},k}$.

Conservative rescaling replaces $q_{m,f}^{\mathrm{rec},k}$ by $\widehat q_{m,f}^k$ in the material correction or its implicit coefficient. It also bounds $\delta F_{z,ms,f}^k$, as derived in Conservative flux rescaling. Without rescaling, $\widehat q_{m,f}^k=q_{m,f}^{\mathrm{rec},k}$ and the formulas above apply directly. With rescaling, their equivalence at convergence holds for the same selected face density and composition correction.

At an absent donor, $q_{m,D}^k=\widehat q_{m,f}^k=0$. Both formulations retain $\dot{\Omega}_f^k z_{ms,D}^{*,k}$ without evaluating a composition or dividing by zero. Material supplied through another face can still pass through the implicit stencil. Including this branch, the implicit material flux used in the species row is

(15.15.1.9)#\[\begin{split}F_{z,ms,f}^{\mathrm{imp},*,k} =\begin{cases} \displaystyle \dot{\Omega}_f^k\frac{\widehat q_{m,f}^k}{q_{m,D}^k} z_{ms,D}^{*,k}+\delta F_{z,ms,f}^k, & q_{m,D}^k>0\\[6pt] \dot{\Omega}_f^k z_{ms,D}^{*,k}, & q_{m,D}^k=0 \end{cases}\end{split}\]

Summing over species cancels the composition correction and gives the material predictor flux:

(15.15.1.10)#\[\begin{split}F_{q,m,f}^{\mathrm{imp},*,k} \equiv\sum_{s\in\mathcal K_m}F_{z,ms,f}^{\mathrm{imp},*,k} =\begin{cases} \displaystyle \dot{\Omega}_f^k\frac{\widehat q_{m,f}^k}{q_{m,D}^k} q_{m,D}^{*,k}, & q_{m,D}^k>0\\[6pt] \dot{\Omega}_f^k q_{m,D}^{*,k}, & q_{m,D}^k=0 \end{cases}\end{split}\]

Here $q_{m,D}^{*,k}=\sum_s z_{ms,D}^{*,k}$. At convergence this becomes $F_{q,m,f}=\dot{\Omega}_f\widehat q_{m,f}$, as does the standard material flux with the same rescaling.

15.15.1.5. Consistent pressure response#

Pressure correction changes the face volume flow. Its material response must follow the chosen species linearization. Hold the predicted densities, the deferred corrections, and the reconstructed ratio fixed while changing $\dot{\Omega}_f$. For a present donor, the responses of the summed species rows are

(15.15.1.11)#\[\left.\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f}\right|_{\mathrm{std}} =q_{m,D}^{*,k}, \qquad \left.\frac{\partial F_{q,m,f}}{\partial\dot{\Omega}_f}\right|_{\mathrm{imp}} =\frac{\widehat q_{m,f}^k}{q_{m,D}^k}q_{m,D}^{*,k}\]

The absent-donor response is $q_{m,D}^{*,k}$ for either formulation. These are responses of the frozen predictor rows, not full derivatives of the converged reconstructed flux. The implicit branch is the response in Eq. (15.7.2.1). The standard branch holds its entire deferred material mass-flow correction fixed along with its composition correction.

15.15.2. Conservative flux rescaling#

Conservative rescaling reduces face corrections before they enter the cell balances. One shared factor is used on both sides of an internal or periodic face, so the flux remains equal and opposite in the neighboring cells. No species inventory is changed independently after the solve.

The basic bound compares a proposed loss with the mass available to support it. For example, if a cell contains $2\ \mathrm{mg}$ of a species at the accepted time level and its outgoing correction would remove $3\ \mathrm{mg}$ during the next step, an inventory-only bound permits at most two-thirds of that correction. Other faces and the implicit solve may supply mass, but this initial bound does not rely on that replenishment.

The construction applies this idea in two stages. First limit the material correction distributed at donor composition. Then limit the composition correction using what remains on the species right-hand side. A bound on their initial sum is insufficient: an outward material correction and an inward composition correction may cancel. Reducing only the latter would increase the net outgoing loss.

15.15.2.1. Common material reconstruction factor#

The material part of Eq. (15.6.5.6) gives the unlimited material correction and its species shares at iterate $k$:

(15.15.2.1)#\[\Delta F_{q,m,f}^k =\dot{\Omega}_f^k\left(q_{m,f}^{\mathrm{rec},k}-q_{m,D_f^k}^k\right), \qquad \Delta F_{z,ms,f}^{q,k} =\Delta F_{q,m,f}^k\widetilde Y_{ms,D_f^k}^{(m),k}\]

Both corrections are zero at an absent donor without evaluating its composition. Their species sum is the material correction. The same inventory calculation is applied to every $q_m$ and every $z_{ms}$. For that calculation, use $c$ to label any one of these components and $Q_{P,c}^n$ for its nonnegative accepted density. Its correction is $\Delta F_q$ for a material or $\Delta F_z^q$ for a species.

For fixed volumes and BDF1, define the total outgoing correction rate $L_{P,c}^k$ and available inventory rate $A_{P,c}^n$:

(15.15.2.2)#\[L_{P,c}^k =\sum_{f\in\mathcal F_P^{\mathrm{int}}} \max\!\left(\sigma_{Pf}\Delta F_{c,f}^k,0\right), \qquad A_{P,c}^n=\frac{\Omega_P}{\Delta t}Q_{P,c}^n\]

The factor $\sigma_{Pf}$ converts the stored face orientation to outward from cell $P$, as in Eq. (15.6.1.3). Taking the positive part counts only correction fluxes that remove mass from $P$. The set $\mathcal F_P^{\mathrm{int}}$ includes internal and periodic faces. A correction can oppose the carrier volume flux, so its own sign determines which cell loses mass.

Both $L$ and $A$ have units of mass per time. In particular, $A$ is the accepted-history contribution from the BDF1 row, rather than the mass itself. The ratio $A/L$ gives the permitted fraction of the outgoing correction:

(15.15.2.3)#\[\begin{split}\eta_{P,c}^k =\begin{cases} 1, & L_{P,c}^k=0\\ \displaystyle\min\!\left(1,\frac{A_{P,c}^n}{L_{P,c}^k}\right), & L_{P,c}^k>0 \end{cases}\end{split}\]

Choose one cell factor that satisfies every component’s bound. Then use the smaller factor from the two cells sharing the face:

(15.15.2.4)#\[\eta_P^k=\min_c\eta_{P,c}^k, \qquad \eta_f^k=\min\!\left(\eta_{\operatorname{cl}(f)}^k, \eta_{\operatorname{cr}(f)}^k\right)\]

Here $0\le\eta_f^k\le1$. A value of one retains the full material correction, and zero removes it. The limited outgoing loss of every component is at most its accepted inventory rate. The same factor applies to every material and species at that face. In floating-point arithmetic, a small relative roundoff margin is reserved in the allowance and active ratios are rounded toward zero. This protects the bound when separate face products are summed; it is not a floor on the transported density.

The resulting material face density is

(15.15.2.5)#\[\widehat q_{m,f}^k =(1-\eta_f^k)q_{m,D_f^k}^k+\eta_f^k q_{m,f}^{\mathrm{rec},k}\]

Thus $\eta_f^k=0$ gives the donor value, while $\eta_f^k=1$ gives the reconstructed value. This nonnegative blend replaces $q^{\mathrm{rec}}$ in the material correction of Eq. (15.6.5.11), or in the coefficient of Eq. (15.15.1.7) (where the reconstructed density is abbreviated as $q_f^k$). Because the factor is the minimum over all components, a species present in only a small amount can restrict the material correction for the entire face.

15.15.2.2. Material-local composition factor#

After fixing the material factor, the proposed composition correction is

(15.15.2.6)#\[\Delta F_{z,ms,f}^{Y,k} =\eta_f^k\dot{\Omega}_f^k q_{m,f}^{\mathrm{rec},k} \delta\widetilde Y_{ms,f}^k\]

Its species sum within each material is zero. The factor $\eta_f^k$ already reduces this correction along with the material reconstruction. The second stage can reduce it further to protect each species.

Recall the species row in Eq. (15.6.7.2), with under-relaxation included. Define ${}^{(z)}\!b_{ms,P,0}^k$ as its assembled right-hand side before adding the composition correction; the subscript $0$ identifies this stage, not a time level. Adding the correction produces

(15.15.2.7)#\[{}^{(z)}\!b_{ms,P}^k ={}^{(z)}\!b_{ms,P,0}^k -\sum_{f\in\mathcal F_P^{\mathrm{int}}} \sigma_{Pf}\delta F_{z,ms,f}^k\]

An outgoing correction is subtracted from the right-hand side; an incoming one adds to it. The remaining allowance is the nonnegative part of ${}^{(z)}\!b_{ms,P,0}^k$ available for outgoing composition correction. It is a mass rate, like the other terms in the species row, and it is not the completed species solution.

This allowance depends on which material linearization was assembled. Let $B_{ms,P}^k$ collect the accepted temporal history, physical sources, boundary inflow, and deferred diffusion terms. Let ${}^{(z)}\!a_{ms,P}^{\mathrm{std},k}$ and ${}^{(z)}\!a_{ms,P}^{\mathrm{imp},k}$ denote the respective unrelaxed diagonals. Using the under-relaxation contribution already described after Eq. (15.6.7.2), with $0<\alpha_z\le1$, gives

(15.15.2.8)#\[\begin{split}\begin{aligned} {}^{(z)}\!b_{ms,P,0}^{\mathrm{std},k} &=B_{ms,P}^k +\frac{1-\alpha_z}{\alpha_z} {}^{(z)}\!a_{ms,P}^{\mathrm{std},k}z_{ms,P}^k -\sum_{f\in\mathcal F_P^{\mathrm{int}}} \sigma_{Pf}\eta_f^k\Delta F_{z,ms,f}^{q,k}\\ {}^{(z)}\!b_{ms,P,0}^{\mathrm{imp},k} &=B_{ms,P}^k +\frac{1-\alpha_z}{\alpha_z} {}^{(z)}\!a_{ms,P}^{\mathrm{imp},k}z_{ms,P}^k \end{aligned}\end{split}\]

The standard row includes the signed material correction. The implicit row already contains it in its matrix. Relaxation must use that same matrix’s diagonal. The allowance therefore accounts for the material, diffusion, and source terms already assembled; the composition correction cannot be assigned the full temporal inventory again. Different linearizations can give different allowances and limiter factors, even though their fluxes agree at convergence when the factors are the same.

Count the proposed outgoing composition corrections in the same way as the material corrections:

(15.15.2.9)#\[L_{P,ms}^{Y,k} =\sum_{f\in\mathcal F_P^{\mathrm{int}}} \max\!\left(\sigma_{Pf}\Delta F_{z,ms,f}^{Y,k},0\right)\]

Within each material, the most restrictive species determines a common composition factor:

(15.15.2.10)#\[\begin{split}\eta_{Y,m,P}^k =\min\!\left( 1,\min_{\substack{s\in\mathcal K_m\\L_{P,ms}^{Y,k}>0}} \frac{\max\!\left({}^{(z)}\!b_{ms,P,0}^k,0\right)}{L_{P,ms}^{Y,k}} \right)\end{split}\]

If there is no positive outgoing composition rate, the factor is one. The face factor and the final composition correction are

(15.15.2.11)#\[\eta_{Y,m,f}^k =\min\!\left(\eta_{Y,m,\operatorname{cl}(f)}^k, \eta_{Y,m,\operatorname{cr}(f)}^k\right)\]

(15.15.2.12)#\[\delta F_{z,ms,f}^k =\eta_{Y,m,f}^k\Delta F_{z,ms,f}^{Y,k}\]

Every species in material $m$ uses the same face factor, preserving the zero sum. Different materials may use different composition factors. The outgoing correction obeys

(15.15.2.13)#\[\sum_{f\in\mathcal F_P^{\mathrm{int}}} \max\!\left(\sigma_{Pf}\delta F_{z,ms,f}^k,0\right) \le\eta_{Y,m,P}^k L_{P,ms}^{Y,k} \le\max\!\left({}^{(z)}\!b_{ms,P,0}^k,0\right)\]

If the pre-correction allowance is nonnegative, the final right-hand side remains nonnegative: the outgoing loss is bounded, and incoming corrections only add to it. The species row explains why this helps. With a positive diagonal, nonnegative neighbor coefficients, and a nonnegative right-hand side, updating a cell from nonnegative neighbor values gives a nonnegative species density. Gauss–Seidel sweeps starting from nonnegative values preserve this property. If those sweeps converge, the linear solution is also nonnegative; the limiter alone does not establish their convergence.

A sufficient bound, not a universal positivity guarantee

A finite negative allowance is not itself a negative species solution. Implicit neighbor inflow can compensate for a negative local right-hand side. The limiter suppresses outgoing composition corrections when that species has no positive allowance, but it cannot guarantee positivity of the whole row in this situation. The completed species predictor still needs an admissibility check. Nonfinite allowances are invalid.

This BDF1 inventory argument does not establish bounds for higher-order time integration, arbitrary sources, or deferred diffusion. It also does not establish nonlinear convergence or the EOS-derived volume-fraction bounds. The face budget above covers internal and periodic corrections. An open boundary with a high-order outgoing correction would need its own matching budget; the displayed bound cannot be extended to it by assumption.

15.15.2.3. Resulting flux and consistency with energy#

At nonlinear convergence, either advection linearization gives

(15.15.2.14)#\[F_{z,ms,f} =\dot{\Omega}_f\widehat q_{m,f}\widetilde Y_{ms,D_f}^{(m)} +\eta_{Y,m,f}\eta_f\dot{\Omega}_f q_{m,f}^{\mathrm{rec}} \delta\widetilde Y_{ms,f}\]

Its material sum is $\dot{\Omega}_f\widehat q_{m,f}$. If $\widehat q_{m,f}>0$, the effective face composition is

(15.15.2.15)#\[\widetilde Y_{ms,f}^{(m)} =\widetilde Y_{ms,D_f}^{(m)} +\frac{\eta_{Y,m,f}\eta_f q_{m,f}^{\mathrm{rec}}}{\widehat q_{m,f}} \delta\widetilde Y_{ms,f}\]

Since $\widehat q_{m,f}\ge\eta_f q_{m,f}^{\mathrm{rec}}$, this is a blend of donor and reconstructed compositions. With both factors one, the unlimited reconstructed flux is recovered. FOU composition and a one-species material have no composition correction.

The mixture face mass flux and advected species enthalpy must use these same limited material and species fluxes, as in Eqs. (15.8.2.2) and (15.8.3.2). Molecular species diffusion remains a separate physical operator. Its enthalpy flux remains the sum of the same physical species diffusion fluxes weighted by species enthalpies, Eq. (15.5.5.4). Multiplying that flux alone by an advection limiter would break the pairing. Any separate diffusion limiter would have to preserve the material’s zero species-flux sum and reuse the resulting species fluxes in energy.

Reducing the correction moves the spatial scheme toward donor-cell transport and can increase interface smearing. The additional cell budgets also require storage and computation. These are costs of enforcing the outgoing-loss bounds, even when the resulting equations converge.

Verification includes conservative, converged transport with FOU species and MIG volume-fraction reconstruction, and with SOU species and MIG where both materials remain present. Sharp-interface tests with SOU species and diffusion, and capillary tests, retain convergence failures. Thus reducing negative-predictor failures does not establish that the coupled nonlinear iteration will converge.

15.15.3. Experimental material-to-energy predictor coupling#

The mixture enthalpy relation in Eq. (15.3.10) is $h=\sum_m w_mh_m$, where $w_m=q_m/\rho$. Even at fixed pressure and temperature, changing the material mass fractions changes $h$. For example, take two materials with enthalpies of $100$ and $300\ \mathrm{kJ/kg}$ at the same pressure and temperature. Increasing the first material’s mass fraction from $0.5$ to $0.6$ changes mixture enthalpy from $200$ to $180\ \mathrm{kJ/kg}$ without changing temperature.

The species predictor changes these mass fractions. The temperature closure then uses the predicted material masses together with the predicted total enthalpy. If the energy predictor has not supplied the corresponding composition-induced enthalpy change, the closure can respond by changing temperature. Consistent face fluxes recover the correct conservative balance at convergence, but independent predictors can still approach it poorly. The experimental update explicitly accounts for the material predictor’s enthalpy change within the energy solve.

Prototype assumptions

The following update has been exercised on a fixed mesh with BDF1 and one species per material. Flow variables and material properties use FOU reconstruction; material volume fractions use FOU or MIG. Species diffusion, reaction, phase change, and conservative flux rescaling are excluded. The energy equation uses the continuity-subtracted form described below. This update is separate from the standard independent predictor sequence.

15.15.3.1. Enthalpy change at fixed pressure and temperature#

First solve the species equations and derive $q_m^{*,k}=\sum_s z_{ms}^{*,k}$. To isolate the effect of the changed material proportions, hold pressure, temperature, and velocity at their incoming values. Evaluate $h_{m,P}^k$ at $(p_P^k,T_P^k)$ and use these same material enthalpies in both mixture sums:

(15.15.3.1)#\[\Delta h_{c,P}^k =\frac{\sum_m q_{m,P}^{*,k}h_{m,P}^k}{\sum_m q_{m,P}^{*,k}} -\frac{\sum_m q_{m,P}^{k}h_{m,P}^k}{\sum_m q_{m,P}^{k}}\]

The subscript $c$ denotes the change due to composition. Here composition means the proportions of materials; each material has only one species in this prototype. Both mixture densities must be positive. Because velocity is held fixed, the kinetic energy $K=\tfrac12\lVert\vec u\rVert^2$ cancels from the difference. Thus $\Delta h_c$ is also the composition contribution to total enthalpy $h_t=h+K$.

The sums require enthalpies only for materials present in the incoming or predicted state. A newly appearing one-species material has a unique composition, so its enthalpy can be evaluated at the incoming pressure and temperature without dividing by its zero incoming mass.

15.15.3.2. The ordinary energy increment equation#

The correction acts on the same assembled energy row used by the ordinary predictor. Using the diagonal and neighbor coefficients of Algebraic Form of the Cell Equation, the action of $A_h^k$ on a cell field $\phi$ is

(15.15.3.2)#\[\left(A_h^k\phi\right)_P ={}^{(h_t)}\!a_P^{\mathrm{rel},k}\phi_P -\sum_{N\in\mathcal N(P)}{}^{(h_t)}\!a_{P,N}^{k}\phi_N\]

The coefficients are those of the energy equation. The diagonal marked $\mathrm{rel}$ includes under-relaxation; the diagonal and neighbor coefficients include the selected implicit transport and numerical enthalpy diffusion. Their units are mass per time. Multiplying by a specific enthalpy therefore gives a cell energy rate. Numerical enthalpy diffusion here is an implicit coupling between neighboring enthalpy increments; its contribution vanishes when those increments vanish. The physical heat-conduction flux remains part of the energy balance.

Write the energy unknown as its incoming value plus an increment, $\Delta h_t=h_t^{*,k}-h_t^k$. Subtracting the matrix action on $h_t^k$ from the predictor row leaves an equation for $\Delta h_t$. Its right-hand side is the negative of the current equation imbalance.

To express that imbalance, evaluate the material balance in Eq. (15.6.1.5) at the incoming iterate. With no material sources, its left-minus-right residual is

(15.15.3.3)#\[R_{m,P}^k =\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\]

It has units of mass per time and is zero when this discrete material balance is satisfied. Its material sum is the mixture continuity residual. Similarly, let $R_{E,P}^k$ be the left-minus-right residual of the conservative energy balance in Eq. (15.5.5.2). It includes pressure time dependence, boundary fluxes, heat transport, and mechanical work, and has units of power. Neither residual includes artificial terms introduced solely to relax the linear solve.

The continuity-subtracted form follows the operation described in Chapter 4: subtract $h_{t,P}^k$ times the discrete continuity residual from the energy residual. This changes the iteration but leaves the converged energy balance unchanged when continuity is satisfied. The ordinary increment equation is therefore

(15.15.3.4)#\[\left(A_h^k\Delta h_t\right)_P=b_{h,P}^k, \qquad b_{h,P}^k=-\left(R_{E,P}^k-h_{t,P}^k\sum_mR_{m,P}^k\right)\]

Here $b_h$ is the right-hand side for the increment, rather than the right-hand side of the equation for the full value $h_t^{*,k}$.

15.15.3.3. Account for the material change in the energy row#

The prototype makes two changes to the ordinary increment equation. First, the completed material predictor supplies the composition contribution $\Delta h_c^k$ from Eq. (15.15.3.1). The energy solve then applies its matrix to the remaining increment $\Delta h_t-\Delta h_c^k$. This remaining part can include pressure and kinetic-energy changes as well as temperature changes.

Second, the ordinary continuity subtraction assigns the same local mixture enthalpy $h_t$ to every material mass imbalance. A mass imbalance in material $m$ instead carries that material’s total enthalpy $h_m+K$. Weighting the material residuals individually removes their associated energy imbalances. Its relation to the ordinary continuity subtraction is the identity

(15.15.3.5)#\[\begin{split}\begin{aligned} R_{E,P}^k-\sum_m(h_{m,P}^k+K_P^k)R_{m,P}^k ={}&\left(R_{E,P}^k-h_{t,P}^k\sum_mR_{m,P}^k\right)\\ &-\sum_m(h_{m,P}^k+K_P^k-h_{t,P}^k)R_{m,P}^k \end{aligned}\end{split}\]

This is a weighted combination of the existing energy and material balances. It has the same converged energy balance when every $R_m$ vanishes. Using this residual to drive the remaining increment is the experimental iteration choice:

(15.15.3.6)#\[\left[A_h^k\left(\Delta h_t-\Delta h_c^k\right)\right]_P =-\left[R_{E,P}^k-\sum_m\left(h_{m,P}^k+K_P^k\right)R_{m,P}^k\right]\]

Expand the matrix action on the left, and move $-A_h^k\Delta h_c^k$ to the right. The residual identity above and Eq. (15.15.3.4) then give $A_h^k\Delta h_t=b_h^k+C_h^k$, where

(15.15.3.7)#\[C_{h,P}^k =\left(A_h^k\Delta h_c^k\right)_P +\sum_m\left(h_{m,P}^k+K_P^k-h_{t,P}^k\right)R_{m,P}^k\]

The term $A_h^k\Delta h_c^k$ applies the same diagonal and neighbor coefficients to the composition increment throughout the mesh. It accounts for both relaxation and numerical transport of that increment. The residual term replaces the mixture-enthalpy weighting of continuity by the individual material enthalpies. Both terms are needed for the displayed equation.

This update creates a real dependency: the material predictor must finish before the energy correction can be evaluated. Thermodynamic properties, transport coefficients, and reconstruction remain frozen at the incoming iterate. It is a segregated Picard iteration with a material solve followed by an energy solve that uses its result.

At nonlinear convergence, $q_m^{*,k}=q_m^k$ and $R_m^k=0$. Then $\Delta h_c^k=C_h^k=0$, and the original conservative energy equation is recovered. At finite tolerance, conservation still depends on the original material and energy residuals. A small residual of the modified linear solve only shows that its chosen increment has been solved accurately. Nonlinear acceptance must still use the ordinary energy measure $b_h$, obtained by removing $C_h$ from the modified right-hand side. The conservative energy residual is related to that measure by

(15.15.3.8)#\[R_{E,P}^k =-b_{h,P}^k+h_{t,P}^k\sum_mR_{m,P}^k\]

Consequently both the unmodified energy measure and the material residuals must be checked.

15.15.3.4. Enthalpy references and scope of the evidence#

An enthalpy reference sets the additive constant in a material’s enthalpy; changing it does not change that material’s heat capacity. For independently conserved materials, consider constant shifts $h_m\mapsto h_m+c_m$. As required by the reference-state identities earlier in this chapter, transform the energy state, histories, boundary fluxes, and sources consistently. Each added $c_m$ multiplies that material’s mass balance, so the conservative energy residual transforms as

(15.15.3.9)#\[R_{E,P}^k\mapsto R_{E,P}^k+\sum_m c_m R_{m,P}^k\]

The sum subtracted from $R_E$ in Eq. (15.15.3.6) gains the same $\sum_m c_m R_m$ term, so the right-hand side is unchanged. The composition increment carries the corresponding change of reference mixture enthalpy. This explains why the update can reduce temperature sensitivity to reference offsets during iteration. Finite linear tolerances and EOS inversion still limit numerical invariance. Independent offsets are not freely interchangeable when mass transfers between materials; reaction and phase-change energy accounting require a separate consistent derivation.

Controlled capillary startup and restart tests show improved convergence under the prototype assumptions. They do not establish curvature accuracy or convergence over the full capillary transient. In a separate translating interface test with a large material-density contrast, thermal equilibrium improved but the energy residual stalled above its acceptance tolerance.

For multispecies materials, the species predictor can change both the material proportions and the composition inside each material. An extension must account for both changes using the material thermodynamic relations, and include every supported source and boundary contribution. Molecular diffusion must retain its species and enthalpy pairing. Conservation, reference-shift behavior, and the original nonlinear residuals remain acceptance requirements for that extension.

15.15.4. Consistent material fractions for MIG reconstruction#

The material fractions used by MIG come from the transported material masses and their EOS densities. At an incoming outer iterate $k$,

(15.15.4.1)#\[\alpha_m^k=\frac{q_m^k}{\rho_m^k} \qquad \alpha_\Sigma^k=\sum_m\alpha_m^k\]

The pressure/material/energy iteration drives $\alpha_\Sigma^k$ toward one. Before convergence, however, the masses and EOS densities can leave a small volume mismatch. MIG receives these intermediate fractions and their gradients when it constructs the next face flux.

The coupled projection in Eq. (15.6.2.14) makes the reconstruction increments sum to zero. If MIG were applied directly to the original cell fractions, adding those increments to the donor fractions would give

(15.15.4.2)#\[\sum_m\alpha_{m,f}^k =\sum_m\left(\alpha_{m,D}^k+\delta\alpha_{m,f}^k\right) =\alpha_{\Sigma,D}^k\]

That would preserve the donor sum, including its current error. The projection also couples the individual material increments. For example, a volume sum error of $6\times10^{-6}$ is small compared with a whole cell, but is about nine percent of a trace fraction of $6.6\times10^{-5}$. A reconstruction change of this scale can substantially alter trace-material transport and feed back into the thermal and pressure updates.

The current MIG scheme instead uses the following temporary fractions that form a consistent partition of the cell. The conserved masses still determine the physical cell fractions through Eq. (15.15.4.1). All quantities in the next formulas use the incoming iterate; superscript $k$ is suppressed.

15.15.4.1. Proportional normalization and its gradient#

Assume $\alpha_m\ge0$ for every material and a finite positive $\alpha_\Sigma$. Define the temporary reconstruction fraction

(15.15.4.3)#\[\beta_m=\frac{\alpha_m}{\alpha_\Sigma} \qquad \sum_m\beta_m=1\]

This preserves relative material proportions and exact absence. Apply this normalization in every cell and at every boundary sample used by the gradient stencil. Then compute the gradient from those normalized values.

For example, a linear discrete gradient in cell $P$ can be written using neighbor-cell differences and boundary-value differences:

(15.15.4.4)#\[(\nabla\beta_m)_P =\sum_{N\in\mathcal N_P}\vec w_{PN}(\beta_{m,N}-\beta_{m,P}) +\sum_{b\in\mathcal B_P}\vec w_{Pb}(\beta_{m,b}-\beta_{m,P})\]

Here $\mathcal N_P$ contains the neighboring cells, $\mathcal B_P$ contains the boundary samples, and the vector weights $\vec w$ depend on the mesh and the chosen gradient stencil. These are the same weights used for other cell gradients. Each normalized sample sums to one across materials, so summing the gradient expression gives

(15.15.4.5)#\[\sum_m(\nabla\beta_m)_P =\sum_{N\in\mathcal N_P}\vec w_{PN}(1-1) +\sum_{b\in\mathcal B_P}\vec w_{Pb}(1-1) =\vec0\]

The identity holds to roundoff. Periodic samples use the corresponding normalized cell values and the usual coordinate transformation.

The continuous quotient rule would instead suggest transforming an existing $\nabla\alpha_m$ locally:

(15.15.4.6)#\[\nabla\!\left(\frac{\alpha_m}{\alpha_\Sigma}\right) =\frac{\nabla\alpha_m-\beta_m\nabla\alpha_\Sigma}{\alpha_\Sigma}\]

Discrete gradients do not generally obey the product or quotient rule. Substituting discrete gradients on the right can therefore give a different result from differentiating the normalized samples directly. Even a small difference can matter when MIG extrapolates a trace fraction to its virtual upstream point. The current scheme uses the direct discrete gradient in Eq. (15.15.4.4).

Use $\beta_m$ and $\nabla\beta_m$ in both adjoining cells throughout the MIG construction: the virtual upstream value, alignment factor, normalized face mapping, and final coupled projection. The upstream bound and donor fallback remain in place. Normalization across materials here precedes the spatial NVD normalization in Eq. (15.6.2.6).

The final projection now acts on increments from $\beta_{m,D}$:

(15.15.4.7)#\[\alpha_{m,f}=\beta_{m,D}+\delta\beta_{m,f} \qquad \sum_m\delta\beta_{m,f}=0 \qquad \sum_m\alpha_{m,f}=1\]

The selected $\alpha_{m,f}$ enters the material face density in Eq. (15.6.3.1) and hence the shared species fluxes. The donor’s original volume sum is not multiplied back into this face result. When $\alpha_\Sigma=1$ throughout the gradient stencil, including its boundary samples, the temporary inputs equal the original ones. Away from that state, normalization changes the reconstructed flux.

15.15.4.2. Complementary treatment of the largest material#

An experimental alternative, not selected in the current scheme, keeps the smaller material fractions and assigns the volume sum correction to the largest one. In each cell, choose the unique largest component $d$. Define

(15.15.4.8)#\[\beta_m=\alpha_m\quad(m\ne d) \qquad \beta_d=1-\sum_{m\ne d}\alpha_m \qquad \nabla\beta_d=-\sum_{m\ne d}\nabla\alpha_m\]

The gradients of the other components remain unchanged. This construction retains their resolved trace variations even when the largest fraction rounds to one. The same principle supplies complementary species differences and gradients in the material diffusion calculation. If the largest component is tied or its complement would be negative, proportional normalization supplies a bounded, material-order-independent fallback for the temporary fractions.

The choice of largest material introduces a switching surface. Consider two almost equal fractions whose sum is $1+\epsilon$. When their ordering reverses, the entire correction $-\epsilon$ moves from one material to the other. The temporary fractions can therefore jump by approximately $|\epsilon|$ even for an arbitrarily small change in the input ordering. An exact-tie fallback does not remove that neighboring jump. Proportional normalization varies smoothly with the fractions while their sum remains positive; the subsequent MIG limiter still has its own branch transitions.

15.15.4.3. Conservation and scope#

Note

These temporary fractions define the face reconstruction. They do not replace the stored species masses, their material sums, or the EOS-derived cell fractions. Other uses of the cell fractions and their gradients, including interface geometry and phase-change models, retain their original inputs. The original cell volume and energy closure residuals must still converge.

Each reconstructed internal-face flux is shared by the adjoining cells with opposite signs, as in Eq. (15.6.1.6). This preserves the finite-volume cancellation of internal fluxes. No species equation is removed, including when a complementary fraction is used in reconstruction.

The procedure addresses consistency of the reconstruction inputs. Conservative flux rescaling instead limits correction losses against a cell’s available mass over the timestep. Neither a unit face-fraction sum nor a bounded virtual upstream fraction alone proves nonnegative solved species or convergence of the coupled iteration.

Controlled slosh continuations show improved convergence with both temporary-fraction constructions. Proportional normalization is the selected method: it avoids dominant-material switching, and an enthalpy-translation comparison converged more reliably with it. Wider tests still require explicit qualification of trace-species cleanup, finite-tolerance conservation, and behavior across materials and meshes.

15.16. References#

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.
J. U. Brackbill, D. B. Kothe, and C. Zemach, “A continuum method for modeling surface tension,” J. Comput. Phys. 100, 335-354 (1992).
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).
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).
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.
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).
C. Kunkelmann, Numerical Modeling and Investigation of Boiling Phenomena, doctoral dissertation, Technische Universität Darmstadt (2011), doi:10.26083/tuprints-00002731.
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.
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).
R. W. Schrage, A Theoretical Study of Interphase Mass Transfer, Columbia University Press (1953).