--- title: "Pressure-Temperature-Equilibrium Multi-Material Algebraic-VOF Framework" --- # Pressure-Temperature-Equilibrium Multi-Material Algebraic-VOF Framework The pressure-temperature-equilibrium (PTE) model describes compressible flow with immiscible materials. All materials present in a cell share one pressure, temperature, and velocity. Each material has its own equation of state (EOS) and composition. Species mix and diffuse within a material. Finite-rate chemistry and transfer between materials enter through conservative sources. Common pressure and temperature do not imply chemical or phase equilibrium. The model advances every species density, mixture momentum, and mixture energy, using total enthalpy as the energy-equation unknown. Material and mixture mass balances follow by summing the species balances. Material volume fractions are obtained from the transported masses and EOS densities. Pressure correction enforces the requirement that the material volumes fill the cell. The finite-volume discretization uses a fixed mesh, the first-order backward differentiation formula (BDF1), algebraic volume-of-fluid (VOF) reconstruction, and a segregated pressure-based iteration. (chap13_materials_and_species)= ## Materials and Species Let $m=1,\ldots,N_m$ index materials and $s=1,\ldots,N_s$ index species. The set $\mathcal K_m$ contains the species indices belonging to material $m$. These sets are disjoint and together contain all species indices. The same chemical constituent can be represented separately in two materials, such as a liquid and its vapor. For a cell of volume $\Omega$, the volume fraction $\alpha_m$ is the fraction occupied by material $m$. An admissible state satisfies ```{math} :label: chap13_alpha_constraints 0\le\alpha_m\le1 \qquad \sum_{m=1}^{N_m}\alpha_m=1 ``` The transported species density $z_{ms}$ is the mass of species $s$ in material $m$ per unit cell volume. Its sum within a material gives the material mass density ```{math} :label: chap13_material_species_mass_sum q_m=\sum_{s\in\mathcal K_m}z_{ms} ``` The mixture density is $\rho=\sum_m q_m$. Both $z_{ms}$ and $q_m$ are mass per unit mixture volume. For a present material, the EOS gives its intrinsic density $\rho_m$, measured per unit volume occupied by that material. Thus ```{math} :label: chap13_cell_masses q_m=\alpha_m\rho_m ``` The volume fraction is recovered as $\alpha_m=q_m/\rho_m$. The species mass fraction within a present material is ```{math} :label: chap13_material_query_composition \widetilde Y_s^{(m)}=\frac{z_{ms}}{q_m} \qquad q_m>0 ``` The vector $\widetilde{\vec Y}^{(m)}$ is the composition supplied to that material's EOS. Its components satisfy ```{math} :label: chap13_material_internal_species \widetilde Y_s^{(m)}\ge0 \qquad \sum_{s\in\mathcal K_m}\widetilde Y_s^{(m)}=1 ``` For a one-species material, $\widetilde Y_s^{(m)}=1$ and $z_{ms}=q_m$. Its species balance is still solved. When a material is absent, $q_m=0$, all its species densities and its volume fraction are zero, and its composition and EOS state are undefined. Mixture sums involving material properties include only present materials. The mass fraction of material $m$ in the mixture is ```{math} :label: chap13_material_mass_fraction w_m=\frac{q_m}{\rho} \qquad \sum_{m=1}^{N_m}w_m=1 ``` For $s\in\mathcal K_m$ in a present material, the mixture species mass fraction is ```{math} :label: chap13_global_from_internal Y_s=\frac{z_{ms}}{\rho}=w_m\widetilde Y_s^{(m)} ``` An absent material has $w_m=0$ and $Y_s=0$ for its species. The mixture species fractions are nonnegative and satisfy ```{math} :label: chap13_vof_species_manifold \sum_{s\in\mathcal K_m}Y_s=w_m \qquad \sum_{s=1}^{N_s}Y_s=1 ``` These fractions describe species amounts in the whole mixture. Material thermodynamics and species diffusion use the composition within each material. ## Solved and Derived Variables The conserved quantities and the unknowns used in the segregated solves are listed below. Here $E$ and $h_t$ denote mixture specific total energy and total enthalpy, respectively. | Equation | Conserved quantity | Solved variable | | --- | --- | --- | | Material species | $z_{ms}$ | Species density $z_{ms}$ | | Momentum | $\rho\vec u$ | Velocity $\vec u$ | | Total energy | $\rho E=\rho h_t-p$ | Total enthalpy $h_t$ | Species sums give the material mass densities $q_m$ and mixture density $\rho$. Composition and mass fractions follow from these densities as defined in {ref}`chap13_materials_and_species`. At the common pressure and temperature, each present material's EOS supplies its intrinsic density $\rho_m$, enthalpy, and other thermodynamic properties. Its volume fraction is then $\alpha_m=q_m/\rho_m$. At fixed pressure, temperature is obtained by matching the mixture enthalpy to the solved total enthalpy minus kinetic energy. Pressure correction enforces that the material volumes fill the cell and corrects the velocity and face volume flux. Transport, pressure, and temperature are coupled by the iteration in {ref}`chap13_segregated_nonlinear_algorithm`. (chap13_material_absence_and_appearance)= ### Material appearance Transport or a physical source can introduce material into a cell where it was absent. The resulting species densities determine its mass and composition, and its EOS density determines its volume fraction. A source must supply consistent species, momentum, and energy contributions; its material mass source is the sum of its species sources, as required by Eq. {eq}`chap13_material_transfer_source_constraint`. Face transport into cells with absent materials is described in {ref}`chap13_material_appearance_face_reference`. ## Thermodynamic Definitions Material thermodynamic properties are evaluated from pressure, temperature, and material composition. The materials share the same cell pressure and temperature, but each material EOS receives only its own composition. It is useful to distinguish a material constitutive model from an evaluated material state. Let $\mathcal{E}_m$ denote the EOS mapping for material $m$. A cell-local thermodynamic state is the result of evaluating that mapping: ```{math} :label: chap13_eos_evaluator_state \mathbf{s}_m(\mathbf{x},t) = \mathcal{E}_m\!\left( p(\mathbf{x},t), T(\mathbf{x},t), \widetilde{\vec{Y}}^{(m)}(\mathbf{x},t) \right). ``` For every active material, the state $\mathbf{s}_m$ supplies finite $\rho_m>0$, $h_m$, $e_m$, the partial specific enthalpies required by species diffusion and phase transfer, the physical sound speed $c_m$ at fixed material composition, and the thermodynamic derivatives required by the PTE closure. It also reports whether the state lies on the material's assigned phase branch and inside its valid thermodynamic domain. The constitutive state is determined only by the local $(p,T,\widetilde{\vec{Y}}^{(m)})$ and the material model $\mathcal E_m$. The evaluated state supplies the intrinsic material density and static thermodynamic enthalpy: ```{math} :label: chap13_material_density \rho_m = \rho_m\!\left( p,T,\widetilde{\vec{Y}}^{(m)};\mathcal{E}_m \right), ``` ```{math} :label: chap13_material_enthalpy h_m = h_m\!\left( p,T,\widetilde{\vec{Y}}^{(m)};\mathcal{E}_m \right). ``` It also supplies the fixed-composition derivatives ```{math} :label: chap13_cmim_material_eos_partial_derivatives \rho_{m,p} = \left(\frac{\partial\rho_m}{\partial p}\right)_{T,\widetilde{\vec{Y}}^{(m)}}, \qquad \rho_{m,T} = \left(\frac{\partial\rho_m}{\partial T}\right)_{p,\widetilde{\vec{Y}}^{(m)}}, \qquad h_{m,p} = \left(\frac{\partial h_m}{\partial p}\right)_{T,\widetilde{\vec{Y}}^{(m)}}, \qquad h_{m,T} = \left(\frac{\partial h_m}{\partial T}\right)_{p,\widetilde{\vec{Y}}^{(m)}}. ``` The mapping $\mathcal E_m$ returns these derivatives with $\rho_m$ and $h_m$, evaluated at the same $(p,T,\widetilde{\vec{Y}}^{(m)})$ on the assigned EOS branch. For the species-diffusion energy flux, let $M_m=q_m\Omega$ be the material mass in the cell and $M_s=z_{ms}\Omega$ the mass of species $s\in\mathcal K_m$. The extensive material enthalpy is $H_m=M_mh_m$. The same evaluated state supplies the partial specific enthalpy ```{math} :label: chap13_species_eos \bar h_{ms} = \left( \frac{\partial H_m}{\partial M_s} \right)_{p,T,M_{\ell\ne s},\,\ell\in\mathcal K_m} \qquad s\in\mathcal{K}_m. ``` Extensivity gives ```{math} :label: chap13_partial_enthalpy_extensivity H_m = \sum_{s\in\mathcal K_m}M_s\bar h_{ms}, \qquad h_m = \sum_{s\in\mathcal K_m}\widetilde Y_s^{(m)}\bar h_{ms}. ``` The partial specific enthalpy is evaluated at the complete material composition with the same reference convention as $h_m$. The effect of consistent enthalpy-reference shifts is derived in {ref}`chap13_enthalpy_references`. Each $\mathcal E_m$ defines the mixing rule for the species inside material $m$. There is no cross-material EOS and no second reconstruction of a material state from mixture species fractions. The mixture volume is the sum of the independently evaluated material volumes: ```{math} :label: chap13_geometric_material_volume \Omega=\sum_m\frac{M_m}{\rho_m}, \qquad \alpha_m=\frac{q_m}{\rho_m}, \qquad \sum_m\frac{q_m}{\rho_m}=1. ``` The material internal energy is related to its static enthalpy by ```{math} :label: chap13_material_internal_energy e_m=h_m-\frac{p}{\rho_m}. ``` The mixture level is the level of the common pressure, velocity, and temperature fields. Each material nevertheless retains its own evaluated EOS state. The mixture density is mass per mixture volume, so each material contributes according to the volume fraction it occupies: ```{math} :label: chap13_mixture_density \rho=\sum_{m=1}^{N_m}q_m=\sum_{m=1}^{N_m}\alpha_m\rho_m ``` The mixture static thermodynamic enthalpy is a specific quantity per unit mixture mass, so it is averaged with the material mass fractions: ```{math} :label: chap13_mixture_enthalpy h=\sum_{m=1}^{N_m}\frac{\alpha_m \rho_m}{\rho} h_m = \sum_{m=1}^{N_m} w_m h_m ``` The mixture internal-energy density follows from the same material masses: ```{math} :label: chap13_mixture_internal_energy \rho e = \sum_{m=1}^{N_m}\alpha_m\rho_m e_m = \rho h-p, ``` where the final equality uses $\sum_m\alpha_m=1$. Defining $K=\frac12\lVert\vec{u}\rVert^2$, the mixture total energy is ```{math} :label: chap13_mixture_total_energy E=e+K, \qquad \rho h_t-p=\rho E. ``` The total enthalpy is defined as the static thermodynamic enthalpy plus kinetic energy: ```{math} :label: chap13_total_enthalpy_def h_t=h+\frac{1}{2}\lVert\vec{u}\rVert^2 ``` ### Pressure-temperature compatibility At a fixed transported state, the species densities $z_{ms}$, their material sums $q_m$, velocity $\vec{u}$, and total enthalpy $h_t$ are known. The materials present in the cell form the set ```{math} :label: chap13_active_material_set \mathcal A=\{m:q_m>0\} ``` For each $m\in\mathcal A$, the species densities determine $\widetilde{\vec{Y}}^{(m)}$ through Eq. {eq}`chap13_material_query_composition`. The material EOS branches are also held fixed. The common pressure and temperature must satisfy the volume and enthalpy compatibility conditions below. With the kinetic energy defined above, the mixture density and target static enthalpy are ```{math} \rho=\sum_{m\in\mathcal A}q_m, \qquad h^{\mathrm{tar}}=h_t-K. ``` A trial $(p,T)$ supplies $\rho_m(p,T,\widetilde{\vec{Y}}^{(m)})$ and $h_m(p,T,\widetilde{\vec{Y}}^{(m)})$ from each active material EOS. Because $q_m\Omega$ is the material mass in a cell, $q_m/\rho_m$ is the fraction of the cell volume occupied by that material. Additive volume requires these fractions to sum to one. Equation {eq}`chap13_mixture_enthalpy` also requires $\sum_{m\in\mathcal A}q_mh_m=\rho h^{\mathrm{tar}}$. These two requirements define the closure residuals ```{math} :label: chap13_pte_closure_residuals F_V(p,T) = \sum_{m\in\mathcal A} \frac{q_m}{\rho_m(p,T,\widetilde{\vec{Y}}^{(m)})}-1 = \sum_{m\in\mathcal A}\alpha_m(p,T)-1 =0, \qquad F_H(p,T) = \sum_{m\in\mathcal A}q_m \left[ h_m(p,T,\widetilde{\vec{Y}}^{(m)})-h^{\mathrm{tar}} \right] =0. ``` Here $\alpha_m(p,T)=q_m/\rho_m(p,T,\widetilde{\vec{Y}}^{(m)})$ is a trial EOS-derived volume fraction, not an independent input. The residual $F_V$ is dimensionless. The residual $F_H$ has units $\mathrm{J/m^3}$, equivalently $\mathrm{Pa}$; dividing it by $\rho$ gives the corresponding specific-enthalpy residual in $\mathrm{J/kg}$. These are compatibility conditions on the transported state, not additional conservation equations. The same residuals enforce the total-energy representation. Using $e_m=h_m-p/\rho_m$ gives ```{math} \left(\sum_{m\in\mathcal A}q_m e_m+\rho K\right) -\left(\rho h_t-p\right) = F_H-pF_V. ``` Therefore, $F_V=F_H=0$ makes the total-energy density reconstructed from the material EOS states equal to the transported total-energy density while the material volumes fill the cell. At the coupled root, this enthalpy form is equivalent to the specific-volume and internal-energy PTE closure discussed in References {ref}`[5] ` and {ref}`[6] `. Assume every active EOS branch is continuous over the common admissible interval and satisfies $c_{p,m}>0$ and $\rho_{m,p}>0$ there. Then ```{math} \left(\frac{\partial F_H}{\partial T}\right)_p = \sum_{m\in\mathcal A}q_m c_{p,m}>0, \qquad \left(\frac{\partial F_V}{\partial p}\right)_T = -\sum_{m\in\mathcal A} \frac{q_m\rho_{m,p}}{\rho_m^2} <0. ``` Thus, at fixed pressure, the enthalpy residual has at most one temperature root. A root exists in the selected interval when the active EOS temperature domains overlap and the endpoint residuals bracket zero. Both residuals generally depend on both $p$ and $T$. For a fixed pressure $p$, the common-temperature inversion is the scalar problem ```{math} :label: chap13_temperature_from_enthalpy f_H(\theta) \equiv \frac{1}{\rho} \sum_{m\in\mathcal A}q_m \left[ h_m(p,\theta,\widetilde{\vec{Y}}^{(m)})-h^{\mathrm{tar}} \right], \qquad f_H(T)=0. ``` During this local solve, $p$, $q_m$, $z_{ms}$, $h^{\mathrm{tar}}$, the active set, and the material EOS branches remain fixed. Each active material EOS is evaluated at every trial temperature $\theta$. CMIM uses this inversion both before and after pressure correction. The pre-pressure evaluation supplies the material densities and derivatives in Eq. {eq}`chap13_fixed_q_volume_response`. The pressure equation uses the local derivative of $F_V$ along $F_H=0$ together with the pressure-induced face-volume-flux response; it is not a local two-variable closure solve. After pressure correction, the corrected velocity defines $h^{\mathrm{tar}}$, the scalar inversion supplies the common temperature, and the material EOS states are evaluated at the corrected $(p,T)$. The complete stage order is given in {ref}`chap13_segregated_nonlinear_algorithm`. (chap13_initial_state)= ### Initial state Prescribe the volume fractions, each present material's composition, and the common pressure, temperature, and velocity. Fractions must satisfy Eqs. {eq}`chap13_alpha_constraints` and {eq}`chap13_material_internal_species`. The common pressure and temperature must lie within every present material's EOS domain on its assigned branch. Evaluate the present material densities and enthalpies, then initialize ```{math} :label: chap13_initial_species_and_enthalpy z_{ms}=\alpha_m\rho_m\widetilde Y_s^{(m)} \qquad h_t=\sum_{m\in\mathcal A}\frac{q_m}{\rho}h_m +\frac12\lVert\vec u\rVert^2 ``` Here $q_m$ and $\rho$ are the species sums. For an absent material, set its species densities to zero without evaluating its composition or EOS state; see {ref}`chap13_material_absence_and_appearance`. These constructions give compatible mass, volume, and energy states, with $\rho h_t-p=\rho E$. Thermodynamic compatibility does not imply mechanical balance. A field intended to remain at rest must also balance pressure against gravity and capillary forces; the conditions are given in {ref}`chap13_equilibria`. ## Mixture Viscosity and Heat Conduction For each material with positive partial density $q_m$, its transport model supplies the dynamic viscosity $\mu_m$ and thermal conductivity $\lambda_m$ at the common cell pressure and temperature and the material-local composition. The cell mixture properties are the normalized volume-fraction averages ```{math} \alpha_\Sigma=\sum_m\alpha_m, \qquad \mu_{\mathrm{mix}}=\frac{1}{\alpha_\Sigma}\sum_m\alpha_m\mu_m, \qquad \lambda_{\mathrm{mix}}=\frac{1}{\alpha_\Sigma}\sum_m\alpha_m\lambda_m. ``` An absent material has $\alpha_m=0$ and does not contribute. At closure, $\alpha_\Sigma=1$; the normalized form is retained during the nonlinear iteration when the volume-fraction sum may differ from unity. For turbulent flow, the momentum and thermal operators use the effective mixture coefficients ```{math} \mu_{\mathrm{eff}}=\mu_{\mathrm{mix}}+\mu_t, \qquad \lambda_{\mathrm{eff}} =\lambda_{\mathrm{mix}}+\frac{\mu_t c_p}{Pr_t}, ``` Here $\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. Turbulent species mixing uses the Schmidt-number closure in Eq. {eq}`chap13_material_effective_diffusivity`. For the Menter SST 2003 closure, the face mass flux uses the CMIM face density and incoming volume flux. Its optional buoyancy correction requires a mixture thermal expansion coefficient. At fixed material mass fractions and material compositions, differentiating $1/\rho=\sum_m w_m/\rho_m$ gives ```{math} :label: chap13_sst_mixture_expansion \beta=-\frac{1}{\rho}\left(\frac{\partial\rho}{\partial T}\right)_{p,w,Y} =\sum_m\alpha_m\beta_m \qquad \beta_m=-\frac{1}{\rho_m}\left(\frac{\partial\rho_m}{\partial T}\right)_{p,Y_m} ``` Present materials supply their EOS derivatives; absent materials contribute zero. Positive buoyancy production enters the $k$ source, while negative production is an implicit sink. This closure represents thermal stratification at fixed composition and excludes production from the material-density jump. Internal-face momentum viscosity is the arithmetic average of the effective cell viscosities. Fourier conduction instead uses a distance-weighted harmonic coefficient. For a face $f$ between cells $L$ and $R$, define ```{math} \theta_f=\frac{d_{Lf}}{d_{Lf}+d_{Rf}}, \qquad \lambda_f^H = \left[ \frac{\theta_f}{\lambda_L} +\frac{1-\theta_f}{\lambda_R} \right]^{-1}, ``` where $d_{Lf}$ and $d_{Rf}$ are the cell-center-to-face distances. Let $\mathcal{G}_f$ denote the positive orthogonal area-over-distance factor in the over-relaxed face decomposition from Eq. {eq}`chap3_scalar_diffusion_flux_split_final`. The two-point conductive energy rate, positive from $L$ to $R$, is ```{math} F_{\lambda,f}^{\mathrm{orth}} = \mathcal{G}_f\lambda_f^H(T_L-T_R). ``` The same $\lambda_f^H$ weights the lagged nonorthogonal correction. The common-temperature closure assumes no unresolved thermal contact resistance, so this Fourier connection remains active across a non-phase-changing material interface. On a link crossing a reconstructed phase-change interface, it is replaced by the two one-sided heat fluxes in Eq. {eq}`chap13_phase_change_energy_closure`. ## Governing Equations Each species has a conservation equation. Summing these equations gives material and mixture mass conservation. Momentum and total energy are mixture balances, using the common velocity, pressure, and temperature. Their domain-integrated forms are given in {ref}`chap13_integral_balances`. ### Species transport For species $s\in\mathcal K_m$ in material $m$, ```{math} :label: chap13_species_transport_target \frac{\partial z_{ms}}{\partial t} +\nabla\cdot(z_{ms}\mathbf u) =-\nabla\cdot\vec{\mathcal J}^{(m)}_s +\alpha_m\dot\omega^{(m)}_s +S^{tr}_{ms} ``` The chemical production rate $\dot\omega_s^{(m)}$ is mass produced per unit material volume per unit time. Multiplication by $\alpha_m$ converts it to a source per mixture volume. The transfer source $S^{tr}_{ms}$ already uses mixture volume and includes reactions between materials, phase change, and external injection or removal. Both source terms in the equation have units $\mathrm{kg/(m^3\,s)}$; positive values add species mass. A transfer source can create an absent material and is not multiplied by its volume fraction. The diffusion flux per unit mixture area is ```{math} :label: chap13_mixture_area_diffusion_flux \vec{\mathcal J}^{(m)}_s=\alpha_m\vec J^{(m)}_s ``` Here $\vec J^{(m)}_s$ is the flux per area occupied by material $m$, in $\mathrm{kg/(m^2\,s)}$. Diffusion redistributes species within a material, so its corrected fluxes satisfy ```{math} :label: chap13_diffusion_constraint_internal \sum_{s\in\mathcal K_m}\vec J^{(m)}_s=\vec0 \qquad \sum_{s\in\mathcal K_m}\vec{\mathcal J}^{(m)}_s=\vec0 ``` A one-species material therefore has zero species diffusion. Reactions within a material conserve its mass: ```{math} :label: chap13_intrinsic_chemistry_mass_constraint \sum_{s\in\mathcal K_m}\dot\omega_s^{(m)}=0 ``` Reactions between materials instead contribute to $S^{tr}_{ms}$ and can change individual material masses. Chemical reactions conserve elemental inventories as well as total mass. Define the material and mixture mass sources by ```{math} :label: chap13_material_transfer_source_constraint S_m=\sum_{s\in\mathcal K_m}S^{tr}_{ms} \qquad \dot m_{tot}=\sum_m S_m ``` ### Material and mixture mass balances Summing Eq. {eq}`chap13_species_transport_target` over the species in material $m$ gives ```{math} :label: chap13_material_continuity_from_species \frac{\partial}{\partial t}\sum_{s\in\mathcal K_m}z_{ms} +\nabla\cdot\left(\mathbf u\sum_{s\in\mathcal K_m}z_{ms}\right) =-\nabla\cdot\sum_{s\in\mathcal K_m}\vec{\mathcal J}^{(m)}_s +\alpha_m\sum_{s\in\mathcal K_m}\dot\omega_s^{(m)} +\sum_{s\in\mathcal K_m}S^{tr}_{ms} ``` The diffusion and chemical production sums within each material vanish. Substituting $q_m=\sum_s z_{ms}$ and the definition of $S_m$ gives ```{math} :label: chap13_material_mass_transport \frac{\partial q_m}{\partial t} +\nabla\cdot(q_m\mathbf u)=S_m ``` Summing over materials, with $\rho=\sum_m q_m$, gives mixture continuity: ```{math} :label: chap13_continuity \frac{\partial\rho}{\partial t} +\nabla\cdot(\rho\mathbf u)=\dot m_{tot} ``` These mass balances follow from the species equations; they are not independent transport equations. Conversion whose removal and addition occur at the same location has $\dot m_{tot}=0$ pointwise. Closed transfer with spatially separated sources instead requires $\int_\Omega\dot m_{tot}\,dV=0$. Its local source remains in continuity. Kunkelmann--Hardt--Wondra redistribution is an example. External injection or removal can change the total mixture mass. ### Mixture momentum ```{math} :label: chap13_momentum \frac{\partial(\rho\mathbf u)}{\partial t} +\nabla\cdot(\rho\mathbf u\otimes\mathbf u) =-\nabla p+\nabla\cdot\boldsymbol\tau +\rho\mathbf g+\mathbf f_\sigma+\mathbf S_u ``` Here $\mathbf g$ is gravitational acceleration, $\mathbf f_\sigma$ is the capillary force per unit volume, and $\mathbf S_u$ is any additional conservative momentum source, in $\mathrm{N/m^3}$. Mass injection or removal requires a specified momentum transfer. A common mixture velocity does not by itself determine that source. With the effective viscosity $\mu$ defined above, the stress is ```{math} :label: chap13_viscous_stress \boldsymbol\tau =\mu\left(\nabla\mathbf u+\nabla\mathbf u^T\right) -\frac23\mu\left(\nabla\cdot\mathbf u\right)\mathbf I ``` For an interface between materials $m$ and $n$, let $\mathbf n_{mn}$ point from $n$ into $m$, and let $\kappa_{mn}$ be positive for a convex region of $m$. The constant surface tension $\sigma_{mn}=\sigma_{nm}$ has units $\mathrm{N/m}$, or equivalently $\mathrm{J/m^2}$. With interfacial area per unit volume $\delta_{\Gamma,mn}$, the continuum surface force is ```{math} :label: chap13_surface_tension_force \mathbf f_{\sigma,mn} =\sigma_{mn}\kappa_{mn}\mathbf n_{mn}\delta_{\Gamma,mn} \qquad \mathbf f_\sigma=\sum_{m0$ is common to all species and materials using this closure, so their turbulent diffusivities are equal within a cell. Laminar flow has $\mu_t=0$. A material modeled without turbulent species mixing has $D_{t,m}=0$. The uncorrected flux per area occupied by material $m$ uses its intrinsic density $\rho_m$: ```{math} :label: chap13_material_raw_fick_flux \widehat{\vec{J}}^{(m)}_s = -\rho_mD^{\mathrm{eff}}_{ms}\nabla\widetilde{Y}^{(m)}_s \qquad s\in\mathcal{K}_m ``` The diffusivities must be valid for the material state, and the mass-fraction-gradient approximation must be appropriate to the mixture. The flux has units of $\mathrm{kg/(m^2\,s)}$. Species-dependent molecular diffusivities can give a nonzero sum of raw fluxes. The corrected fluxes satisfy Eq. {eq}`chap13_diffusion_constraint_internal`. The face treatment in Eq. {eq}`chap13_material_corrected_face_flux` enforces this constraint by rescaling positive and negative raw fluxes separately within each material. (chap13_species_discretization)= ## Species Discretization CMIM advances the species densities $z_{ms}$. Their sum within each material is its mass per mixture volume, $q_m=\sum_{s\in\mathcal K_m}z_{ms}$, so species transport must also conserve material mass. For a fixed cell $P$, integration of Eq. {eq}`chap13_species_transport_target` and BDF1 time discretization give ```{math} :label: chap13_species_bdf1_balance \frac{\Omega_P}{\Delta t} \left(z_{ms,P}^{n+1}-z_{ms,P}^{n}\right) +\sum_{f\in\mathcal F_P}\sigma_{Pf} \left(F_{z,ms,f}^{n+1}+F_{ms,f}^{n+1}\right) =\Omega_P S_{z,ms,P}^{n+1} ``` Here $F_{z,ms,f}$ and $F_{ms,f}$ are the advective and diffusive mass flow rates, both in $\mathrm{kg/s}$. Internal-face fluxes are positive from $L=\operatorname{cl}(f)$ to $R=\operatorname{cr}(f)$; the incidence sign $\sigma_{Pf}$ converts them to cell-outward fluxes as defined in {ref}`app_loci_face_orientation`. The source per mixture volume is $S_{z,ms}=\alpha_m\dot\omega_s^{(m)}+S_{ms}^{tr}$, and $n$ denotes the accepted time level. ### Face fluxes The face volume flow rate $\dot\Omega_f=\vec u_f\cdot\vec S_f$ comes from momentum interpolation. It selects the donor $D$: $L$ for positive flow and $R$ for negative flow. For a material present in the donor, reconstruct its volume fraction, intrinsic density, and composition separately, then form ```{math} :label: chap13_reconstructed_species_flux q_{m,f}=\alpha_{m,f}\rho_{m,f} \qquad z_{ms,f}=q_{m,f}\widetilde Y_{s,f}^{(m)} \qquad F_{z,ms,f}=\dot\Omega_f z_{ms,f} ``` Every species of a material uses the same $q_{m,f}$. Its face composition sums to one, so the species fluxes sum to the material mass flow: ```{math} :label: chap13_q_face_flux \sum_{s\in\mathcal K_m}F_{z,ms,f} =F_{q,m,f}=\dot\Omega_f q_{m,f} ``` First-order upwinding (FOU) uses donor values. Second-order upwinding (SOU) extrapolates limited slopes. Modified IG (MIG) sharpens the volume-fraction profile; density and composition still use their selected FOU or SOU reconstruction. A conditional material state uses SOU only when the material is present throughout the donor's gradient stencil. Otherwise it retains the donor value. An absent donor supplies zero reconstructed material and species fluxes without evaluating its composition or EOS state. The formulas are given in {ref}`chap13_face_reconstruction`. Diffusion redistributes species within a material. It requires that material to be present in both adjacent cells and gives zero flux for a one-species material. The corrected Fick flux in {ref}`chap13_species_diffusion_flux` satisfies ```{math} :label: chap13_material_species_diffusion_flux_sum \sum_{s\in\mathcal K_m}F_{ms,f}=0 ``` Thus diffusion contributes no net material mass flow. The same physical species flux enters the diffusive enthalpy transport. Phase change is represented by transfer sources, not by cross-material Fick diffusion. (chap13_species_sources_and_predictor)= ### Implicit species update At nonlinear iterate $k$, solve for the predictor $z_{ms}^{*,k}$ while holding the volume flow, reconstructed face states, material properties, and source coefficients fixed. For advection, keep the donor contribution implicit and evaluate the face-minus-donor correction at the incoming state: ```{math} :label: chap13_standard_species_predictor_flux F_{z,ms,f}^{*,k} =\dot\Omega_f^k z_{ms,D}^{*,k} +\dot\Omega_f^k\left(z_{ms,f}^k-z_{ms,D}^k\right) ``` The correction vanishes for FOU. At nonlinear convergence, $z_{ms,D}^{*,k}=z_{ms,D}^k$, and the expression recovers the reconstructed flux. At an absent donor the known correction is zero, but the implicit donor term remains: material supplied through another face can pass through the stencil. A bounded reconstruction alone does not ensure a nonnegative predictor; the outgoing correction can exceed the mass available in a cell. For diffusion, retain the corrected two-point conductance $G_{ms,f}^k$ implicitly and hold the material densities $q_m^k$ fixed. On a supported link, ```{math} :label: chap13_material_linearized_diffusion_flux F_{ms,f}^{\mathrm{lin},*,k} =F_{ms,f}^k+G_{ms,f}^k \left[ \frac{z_{ms,L}^{*,k}-z_{ms,L}^k}{q_{m,L}^k} -\frac{z_{ms,R}^{*,k}-z_{ms,R}^k}{q_{m,R}^k} \right] ``` The conductance is defined in {ref}`chap13_species_diffusion_flux`. Unsupported links give zero without evaluating either ratio. The frozen $q_m^k$ makes this a linear approximation in species density; the ratio $z_{ms}^{*,k}/q_m^k$ is not the predicted composition. The bracket vanishes at nonlinear convergence, leaving the physical flux $F_{ms,f}$. Sources use their incoming values, with a nonnegative sink coefficient $\ell_{ms}^k$ where a destruction term is treated implicitly: ```{math} :label: chap13_species_linearized_source S_{z,ms}^{\mathrm{lin},*,k} =S_{z,ms}^k-\ell_{ms}^k\left(z_{ms}^{*,k}-z_{ms}^k\right) ``` Explicit sources have $\ell_{ms}^k=0$. Substituting these fluxes and sources into the time-discrete balance gives the species predictor: ```{math} :label: chap13_discrete_direct_species_equation \frac{\Omega_P}{\Delta t} \left(z_{ms,P}^{*,k}-z_{ms,P}^{n}\right) +\sum_{f\in\mathcal F_P}\sigma_{Pf} \left(F_{z,ms,f}^{*,k}+F_{ms,f}^{\mathrm{lin},*,k}\right) =\Omega_P S_{z,ms,P}^{\mathrm{lin},*,k} ``` Each species, including the sole species of a one-species material, has an assembled row of the form ```{math} :label: chap13_species_algebraic_row {}^{(z)}\!a_{ms,P}^{k}z_{ms,P}^{*,k} =\sum_{N\in\mathcal N(P)} {}^{(z)}\!a_{ms,P,N}^{k}z_{ms,N}^{*,k} +{}^{(z)}\!b_{ms,P}^{k} ``` The diagonal contains the time term, outgoing convection, two-point diffusion, and $\Omega_P\ell_{ms,P}^k$. Neighbor coefficients contain incoming convection and two-point diffusion. The right-hand side contains the accepted history, known sources, and deferred flux corrections. For species under-relaxation $0<\alpha_z\le1$, divide the diagonal by $\alpha_z$ and add $(1-\alpha_z)\,{}^{(z)}\!a_{ms,P}^kz_{ms,P}^k/\alpha_z$ to the right-hand side, using the unrelaxed diagonal in this addition. The added term cancels the diagonal change at nonlinear convergence. {ref}`chap13_chemistry_integration` defines the chemistry sink coefficient and the alternative split treatment. Split chemistry replaces the accepted species history by the reaction endpoint and omits those reactions from the transport source. {ref}`chap13_material_advection_variants` moves the material part of the deferred correction into the donor coefficient. {ref}`chap13_conservative_rescaling` bounds the outgoing face corrections. ### Material conservation The material predictor is always the species sum, $q_m^{*,k}=\sum_s z_{ms}^{*,k}$; it requires no separate transport solve. At nonlinear convergence, summing the species equations cancels physical diffusion and reactions within each material, giving ```{math} :label: chap13_q_bdf1_balance \frac{\Omega_P}{\Delta t} \left(q_{m,P}^{n+1}-q_{m,P}^{n}\right) +\sum_{f\in\mathcal F_P}\sigma_{Pf}F_{q,m,f}^{n+1} =\Omega_P S_{m,P}^{n+1} ``` Here $S_m=\sum_s S_{ms}^{tr}$ is the net transfer into the material. Each shared face flux enters its neighboring cells with opposite signs. For split chemistry, the same balance uses the reacted history $q_m^{\mathrm{chem}}=\sum_s z_{ms}^{\mathrm{chem}}$ and only the remaining transfer sources. During the nonlinear iteration, the implicit diffusion and chemistry sink corrections need not sum to zero within a material. Under-relaxation also adds a temporary term. These terms vanish at nonlinear convergence; a converged linear solve alone does not establish the material balance. The convergence checks are given in {ref}`chap13_segregated_nonlinear_algorithm`. The species-derived material masses then supply the volume-closure error used by pressure correction. (chap13_pressure_correction)= ## Pressure-Correction Equation The species predictor determines the material masses, but their EOS-derived volumes may not fill the cell. Pressure correction couples the change in these volumes to the change in face transport. The species densities remain fixed during the pressure solve; subsequent nonlinear iterations bring their transport into agreement with the corrected flow. ### Volume error and thermodynamic response At pressure $p_P^k$, solve Eq. {eq}`chap13_temperature_from_enthalpy` using the predicted species densities, total enthalpy, and velocity. This gives the common temperature $T_P^{\mathrm{pred},k}$ and intrinsic material densities $\rho_{m,P}^{\mathrm{pred},k}$. The predicted volume error is ```{math} :label: chap13_fixed_q_volume_error 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 ``` A positive value means the materials require more than the available cell volume; a negative value means they occupy less. Absent materials contribute zero and require no intrinsic density. Pressure changes both density and the temperature needed to preserve the predicted static enthalpy $h=h_t^{*,k}-\tfrac12|\vec u^{*,k}|^2$. For the following two local derivatives, omit cell and iteration labels and let $q_m$ denote the fixed predictor value. All EOS quantities are evaluated at the pre-pressure state, and sums include only present materials. Differentiating $\sum_m q_mh_m=\rho h$ at fixed species densities and $h$ gives ```{math} :label: chap13_fixed_q_temperature_pressure_response \left(\frac{dT}{dp}\right)_{h,z} =-\frac{\displaystyle\sum_m q_m h_{m,p}} {\displaystyle\sum_m q_m c_{p,m}} ``` Here $c_{p,m}$ is the material heat capacity and $h_{m,p}=1/\rho_m+T\rho_{m,T}/\rho_m^2$. The EOS derivatives $\rho_{m,p}$ and $\rho_{m,T}$ hold composition fixed, and hold temperature or pressure fixed, respectively. The resulting volume response is ```{math} :label: chap13_fixed_q_volume_response \chi =-\left.\frac{dE_V}{dp}\right|_{h,z} =\sum_m\frac{q_m}{\rho_m^2} \left[\rho_{m,p}+\rho_{m,T}\left(\frac{dT}{dp}\right)_{h,z}\right] ``` The linearization requires $\sum_m q_m c_{p,m}>0$ and $\chi>0$. The coefficient $\chi$ has units of $\mathrm{Pa^{-1}}$: a pressure increase $p'$ reduces the local volume error by $\chi p'$ to first order. ### Pressure response of face transport The species predictor uses the incoming volume flow $\dot\Omega_f^k$. During pressure correction, hold its donor $D$, predicted species densities, reconstruction coefficients, and deferred mass-flow corrections fixed. For standard deferred advection, only the implicit donor contribution in Eq. {eq}`chap13_standard_species_predictor_flux` changes with volume flow: ```{math} :label: chap13_fixed_q_face_flux_derivative \left.\frac{\partial F_{q,m,f}}{\partial\dot\Omega_f}\right|_{\mathrm{std}} =q_{m,D}^{*,k} ``` For implicit material advection, multiply this response by $\widehat q_{m,f}^k/q_{m,D}^k$ when the incoming donor contains the material; use the same donor response when it is absent. Here $\widehat q_{m,f}^k$ is the selected material face density, including conservative rescaling when used. Both responses are derived in Eq. {eq}`chap13_advection_variant_pressure_response`. For cell $P$, divide each material's mass-flow response by its predicted intrinsic density and sum: ```{math} :label: chap13_fixed_q_receiver_face_weight W_{P,f}^k =\sum_{\substack{m\\q_{m,P}^{*,k}>0}} \frac{1}{\rho_{m,P}^{\mathrm{pred},k}} \frac{\partial F_{q,m,f}}{\partial\dot\Omega_f} ``` This dimensionless weight converts a change in face volume flow into a change in the material volume flow measured at cell $P$'s densities. The two cells sharing a face can therefore have different weights, although the material mass flux is shared. Their material sets and EOS densities remain fixed in this linearization. Let $L=\operatorname{cl}(f)$ and $R=\operatorname{cr}(f)$. The selected flux formulation supplies the compact pressure correction ```{math} :label: chap13_fixed_q_face_volume_flux_correction \dot\Omega_f^{\prime,k} =d_f^{\mathrm{vol}}\left(p_L^{\prime,k}-p_R^{\prime,k}\right) ``` The coefficient $d_f^{\mathrm{vol}}$, in $\mathrm{m^3/(s\,Pa)}$, includes the formulation's SIMPLEC response and any direct pressure dependence of its numerical mass flux. It is the complete coefficient used in the compact pressure row and in the face-flow update. Its construction depends on the selected formulation; the common cell-momentum approximations are derived in {numref}`Chapter %s `. ### Pressure-correction balance The momentum predictor supplies $\dot\Omega_f^{*,k}$, while the species predictor used $\dot\Omega_f^k$. The total change in transport seen by the pressure equation is their difference plus $\dot\Omega_f^{\prime,k}$. For fixed cell volumes and BDF1, approximate the associated material-volume change by the net additional outward volume flow over $\Delta t$. Using the frozen weights above, the linearized volume balance is ```{math} :label: chap13_linearized_volume_balance E_{V,P}^{\mathrm{pred},k}-\chi_P^k p_P^{\prime,k} -\frac{\Delta t}{\Omega_P} \sum_{f\in\mathcal F_P}\sigma_{Pf}W_{P,f}^k \left(\dot\Omega_f^{*,k}-\dot\Omega_f^k+\dot\Omega_f^{\prime,k}\right) =0 ``` The first correction accounts for local compression at fixed species masses. The second estimates the effect of changed transport: additional outflow reduces the occupied volume. This is a segregated approximation; it does not re-solve the species equations or change their stored predictors. The incidence sign $\sigma_{Pf}$ converts the stored face orientation to cell-outward flow. Multiply by $\Omega_P/\Delta t$ and substitute the face-pressure correction. For an interior cell, let $N_f$ be the neighbor across face $f$. Collecting the unknown pressure corrections gives :::{admonition} CMIM pressure-correction equation ```{math} :label: chap13_fixed_q_pressure_row \begin{aligned} \left(\frac{\Omega_P}{\Delta t}\chi_P^k +\sum_{f\in\mathcal F_P}W_{P,f}^k d_f^{\mathrm{vol}}\right)p_P^{\prime,k} -\sum_{f\in\mathcal F_P}W_{P,f}^k d_f^{\mathrm{vol}}p_{N_f}^{\prime,k} \\ =\frac{\Omega_P}{\Delta t}E_{V,P}^{\mathrm{pred},k} -\sum_{f\in\mathcal F_P}\sigma_{Pf}W_{P,f}^k \left(\dot\Omega_f^{*,k}-\dot\Omega_f^k\right) \end{aligned} ``` The same face response weights the predictor-flow difference and the pressure correction. Physical mass sources already enter through the species predictor and its volume error; adding them again would count them twice. ::: Periodic faces use the same coupling. A prescribed volume flux has no pressure correction; stationary impermeable boundaries have zero normal volume flow. At a prescribed-pressure boundary, use its specified $p'$ in the face correction and the boundary donor state in the material weight. ### Prescribed mass and volume rates Let $A_b=\sum_{f\in b}A_f$ be the area of a fixed boundary patch, and let $\dot m_b^{\mathrm{in}}$ be its prescribed mass rate, positive into the domain. The outward mass flux per unit area and the incoming material state satisfy ```{math} :label: chap13_prescribed_mass_boundary j_b=-\frac{\dot m_b^{\mathrm{in}}}{A_b},\qquad q_{m,f}=\alpha_{m,b}\rho_m(p_f,T_b,Y_{m,b}),\qquad \rho_f=\sum_m q_{m,f},\qquad \dot\Omega_f^*=\frac{j_b A_f}{\rho_f}. ``` The entering species mass rate is $\dot m_{ms,f}=j_b A_f(q_{m,f}/\rho_f)Y_{ms,b}$. Thus material mass fractions, rather than volume fractions, partition the prescribed total mass rate. Exactly absent materials require no EOS query. For outward flow the interior transport donor supplies the density and composition. A prescribed volume rate instead sets $\dot\Omega_f^*=-\dot\Omega_b^{\mathrm{in}}A_f/A_b$ directly. Within a fixed-mass pressure solve, the selected boundary density is frozen. The prescribed-rate boundary therefore has $\dot\Omega_f^{\mathrm c}=\dot\Omega_f^*$ and no boundary pressure conductance. Its predictor difference $\dot\Omega_f^*-\dot\Omega_f^k$ still enters {eq}`chap13_fixed_q_pressure_row`. The next nonlinear iteration refreshes the donor state and density. At convergence, the material/species transport and the velocity boundary carry the same prescribed mass rate. A fixed-mass outlet can retain the interior normal-speed shape $s_f=|\vec u_i\cdot\vec n_f|$. With outward total rate $M_b\ge0$, set ```{math} :label: chap13_fixed_mass_outlet_shape D_b=\sum_{f\in b}\rho_f s_f A_f,\qquad u_{n,f}=\frac{M_b s_f}{D_b},\qquad \sum_{f\in b}\rho_f u_{n,f}A_f=M_b. ``` When $D_b=0$, use $u_{n,f}=M_b/(\rho_f A_b)$. In both cases the tangential interior velocity is preserved. The density in this normalization is the selected transport density, including a conservative high-order correction when that correction is active. ### Corrected state After solving for $p^{\prime,k}$, apply pressure relaxation and bounds: ```{math} :label: chap13_cmim_pressure_under_relaxation p_P^{\mathrm c,k} =\operatorname{clip}\!\left(p_P^k+\alpha_p p_P^{\prime,k};p_{\min},p_{\max}\right) ``` Here $\alpha_p$ is the pressure relaxation factor, and $\mathrm c,k$ denotes the corrected state. The cell velocity correction is ```{math} :label: chap13_cmim_velocity_pressure_update \vec u_P^{\mathrm c,k} =\vec u_P^{*,k} -\frac{D_P}{\Omega_P} \sum_{f\in\mathcal F_P}\sigma_{Pf}p_f^{\prime,k}\vec S_f ``` The cell response $D_P$ includes momentum relaxation and the SIMPLEC approximation. The face value $p_f^{\prime,k}$ uses the pressure-force approximation of the selected formulation. These are the cell corrections from {numref}`Chapter %s `; the face volume flow uses the matching compact response: ```{math} :label: chap13_cmim_volume_flux_pressure_update \dot\Omega_f^{\mathrm c,k} =\dot\Omega_f^{*,k} +d_f^{\mathrm{vol}}\left(p_L^{\prime,k}-p_R^{\prime,k}\right) ``` At the corrected pressure, solve the common-temperature closure again with the unchanged species predictors and static enthalpy $h_t^{*,k}-\tfrac12|\vec u^{\mathrm c,k}|^2$. Reevaluate material densities and volume fractions. A single pressure solve need not close both transport and thermodynamics; the checks in {ref}`chap13_pressure_transport_checks` determine whether another nonlinear iteration is needed. (chap13_consistent_transport)= ## Consistent Transport The species mass fluxes in {ref}`chap13_species_discretization` determine what mass crosses each face. Momentum and energy transport must use these same fluxes, including any limits applied to the reconstruction. The contact states in {ref}`chap13_equilibria` give limiting solutions for this coupled transport. ### Mass and momentum At incoming nonlinear iterate $k$, the mixture face density is $\rho_f^k=\sum_m\widehat q_{m,f}^k$. Here $\widehat q_{m,f}$ is the reconstructed $q_{m,f}$, or its limited value from Eq. {eq}`chap13_separate_material_face_value`. Summing the species fluxes gives the mixture mass flow: ```{math} :label: chap13_incoming_picard_mass_flux \dot m_f^k =\rho_f^k\dot\Omega_f^k =\sum_m F_{q,m,f}^k =\sum_m\sum_{s\in\mathcal K_m}F_{z,ms,f}^k ``` This mass flow carries momentum with the selected FOU or SOU velocity: ```{math} :label: chap13_material_momentum_flux_identity \vec F_{\rho u,f}^{\mathrm{adv},k} =\dot m_f^k\vec u_f^{\mathrm{adv},k} ``` The transported velocity $\vec u_f^{\mathrm{adv},k}$ need not equal the pressure-coupled face velocity that determines $\dot\Omega_f^k$. ### Advected enthalpy Each species mass flux carries its partial specific enthalpy. Let $\bar h_{ms,f}^{\mathrm{adv},k}$ be its FOU donor value or SOU reconstruction from the incoming cell enthalpies. SOU uses the same material-support condition as the species reconstruction; otherwise it retains the donor value. Adding the kinetic energy carried by the mixture mass flow gives ```{math} :label: chap13_incoming_direct_z_enthalpy_flux F_{\rho h_t,f}^{\mathrm{adv},k} =\sum_m\sum_{s\in\mathcal K_m} F_{z,ms,f}^k\bar h_{ms,f}^{\mathrm{adv},k} +\frac12\dot m_f^k\lVert\vec u_f^{\mathrm{adv},k}\rVert^2 ``` A material absent from the donor contributes zero without evaluating its enthalpy. With FOU, the expression uses donor composition, enthalpies, and velocity. With SOU, independently extrapolating mixture enthalpy generally does not give this species-weighted flux. A composition correction can sum to zero mass transport and still carry enthalpy because the species enthalpies differ. Its contribution must therefore remain in the energy flux. Any conservative rescaling must also use the same limited species fluxes in both equations. ### Diffusive enthalpy Species diffusion also carries enthalpy. The energy equation uses the same physical diffusion flux $F_{ms,f}^k$ as the species equation: ```{math} :label: chap13_cmim_species_enthalpy_face_flux F_{hY,f}^k =\sum_m\sum_{s\in\mathcal K_m} \bar h_{ms,f}^{\mathrm{diff},k}F_{ms,f}^k ``` The face enthalpy for diffusion is distance-centered: ```{math} :label: chap13_cmim_species_enthalpy_face_value \bar h_{ms,f}^{\mathrm{diff},k} =(1-\theta_f)\bar h_{ms,L}^k+\theta_f\bar h_{ms,R}^k ``` Here $\theta_f=d_{Lf}/(d_{Lf}+d_{Rf})$ uses the cell-center-to-face distances. This interpolation is independent of the advective upwind direction. Unsupported material links contribute zero, as specified in {ref}`chap13_species_diffusion_flux`. The species diffusion fluxes sum to zero within each material. They carry no net material mass or common kinetic energy, so static species enthalpies suffice. Their enthalpy-weighted sum can still be nonzero. ### Energy update With these fluxes, BDF1 discretization of Eq. {eq}`chap13_cmim_total_enthalpy_balance` on a fixed cell gives ```{math} :label: chap13_total_enthalpy_bdf1_pressure_time \begin{aligned} &\frac{\Omega_P}{\Delta t} \left[(\rho h_t)_P^{n+1}-(\rho h_t)_P^n\right] +\sum_{f\in\mathcal F_P}\sigma_{Pf} \left(F_{\rho h_t,f}^{\mathrm{adv},n+1}+F_{hY,f}^{n+1}\right) \\ &\qquad =\frac{\Omega_P}{\Delta t} \left(p_P^{n+1}-p_P^n\right)+\mathcal R_{E,P} \end{aligned} ``` The incidence sign $\sigma_{Pf}$ gives each internal face opposite signs in its adjacent cells. The remaining energy rate $\mathcal R_{E,P}$ contains conduction, viscous work, gravity and capillary work, and conservative sources, evaluated at the new time level. The pressure-time term uses the incoming pressure during a segregated iteration. For the advective flux, retain donor enthalpy implicitly and hold the remaining contribution at the incoming iterate: ```{math} :label: chap13_material_enthalpy_deferred_correction F_{\rho h_t,f}^{\mathrm{lin},*,k} =\dot m_f^k h_{t,D}^{*,k} +\left(F_{\rho h_t,f}^{\mathrm{adv},k} -\dot m_f^k h_{t,D}^k\right) ``` The bracket is fixed during the linear solve. At nonlinear convergence, $h_{t,D}^{*,k}=h_{t,D}^k$, leaving the required species-weighted flux. The species, momentum, and energy equations share the incoming face reconstruction, but their different predictor unknowns need not give identical mass transport within an unfinished nonlinear iteration. The energy update also accounts for the completed species predictor's composition change at fixed incoming pressure and temperature. Together with an enthalpy-weighted species-residual correction, this couples the species and energy updates without changing the converged conservation balance. The derivation is given in {ref}`chap13_material_energy_coupling`. (chap13_surface_tension_and_curvature)= ## Surface Tension and Curvature Surface tension supplies the capillary force in Eq. {eq}`chap13_momentum` and its work in Eq. {eq}`chap13_cmim_total_enthalpy_balance`. Curvature is estimated from the material fractions; it is not a transported bulk variable. The force model assumes constant surface tension for each immiscible material pair. It therefore excludes tangential stresses caused by variations in surface tension, known as Marangoni stresses. ### Capillary force For two materials $m,n$, define a color from their volume fractions: ```{math} :label: chap13_celeste_pair_color C=C_{mn}=\frac{\alpha_m}{\alpha_m+\alpha_n} \qquad \alpha_m+\alpha_n>0 ``` Color one identifies material $m$ and color zero identifies material $n$. At volume closure in a two-material cell, $C=\alpha_m$. Write $\sigma=\sigma_{mn}$ for the pair tension. The normal points toward increasing color, and curvature is its negative divergence: ```{math} :label: chap13_celeste_curvature_fit \vec n=\frac{\nabla C}{|\nabla C|} \qquad \kappa=-\nabla\cdot\vec n \qquad |\nabla C|>0 ``` For a convex drop of material $m$, the normal points inward and curvature is positive. The continuum-surface-force model {ref}`[2] ` spreads the interfacial force over the cells where color changes. With $\vec n\delta_\Gamma=\nabla C$, the force density in Eq. {eq}`chap13_surface_tension_force` becomes ```{math} :label: chap13_pair_capillary_force \vec f_\sigma=\sigma\kappa\nabla C ``` Here $\sigma$ has units of $\mathrm{N/m}$, curvature has units of inverse length, and $\vec f_\sigma$ has units of $\mathrm{N/m^3}$. Reversing the material pair reverses both curvature and the color gradient, leaving the force unchanged. For a circular cylinder or sphere of radius $R$, the Young--Laplace equilibrium is ```{math} :label: chap13_surface_tension_laplace_jump \kappa=\frac{d-1}{R} \qquad p_m-p_n=\sigma\kappa ``` Here $d=2$ for a circular cylinder and $d=3$ for a sphere. A junction of three materials requires the extension in {ref}`chap13_three_material_interfaces`. ### Curvature methods All curvature methods supply $\kappa$ to the same force calculation. They differ in how they estimate geometry from cell-average fractions. | Method | Curvature estimate | Main approximation | | --- | --- | --- | | {ref}`Direct divergence ` | Differentiate color, interpolate gradients to faces, and sum normal fluxes | Weak gradients require regularization | | {ref}`Smoothed face or nodal divergence ` | Smooth a geometry color, form face or vertex normals, and integrate their divergence | Smoothing width and normal averaging depend on the mesh | | {ref}`CELESTE ` | Fit color derivatives and normal derivatives, then average curvature | Local Taylor fits require adequate neighboring data | | {ref}`Geometric fit ` | Reconstruct interface planes and fit a paraboloid to them | Requires convex cells with planar faces and enough interface pieces for the fit | | {ref}`Auxiliary level set (RLSG) ` | Reconstruct a distance field on a Cartesian grid and transfer its curvature to the flow cells | Both interface reconstruction and auxiliary-grid spacing affect the result | | Prescribed curvature | Supply a constant value | Does not follow changes in interface shape | Smoothing acts on a separate geometry field. The force still uses the original color to locate the interface, and the transported fractions remain unchanged. Regularization suppresses normals where the color gradient is small; it can also reduce the estimated curvature. Neither operation ensures curvature accuracy under mesh refinement. Geometric fits and RLSG likewise leave material transport to the species equations. The derivations in {ref}`chap13_curvature_derivations` use fixed meshes in Cartesian coordinates and give the geometry and wall assumptions for each method. (chap13_capillary_balance)= ### Pressure balance and work In the absence of other forces, a stationary interface requires pressure to balance the capillary force. A balanced-force discretization {ref}`[3] ` uses compatible pressure and color operators. The cell force is ```{math} :label: chap13_capillary_cell_force \vec f_{\sigma,P} =\frac{\sigma\kappa_P}{\Omega_P} \sum_{f\in\mathcal F_P}\sigma_{Pf}C_f\vec S_f ``` The face color uses the interpolation weights of the selected pressure formulation. For constant curvature and a matching pressure reconstruction, $p_P=p_0+\sigma\kappa C_P$ gives $p_f=p_0+\sigma\kappa C_f$. Substitution into the pressure-force sum then cancels the capillary force on a closed cell, provided the boundary values satisfy the same equilibrium. The face-velocity calculation must preserve this balance too. With no other forces, its pressure difference is balanced when ```{math} :label: chap13_capillary_face_balance p_R-p_L-\sigma\kappa_f(C_R-C_L)=0 ``` Here $L=\operatorname{cl}(f)$ and $R=\operatorname{cr}(f)$. Curvature is interpolated with the matching face weights. Momentum interpolation accounts for this capillary pressure change and removes the interpolated cell-force contribution already present in the velocity predictor. Two errors can drive spurious currents: a mismatch between these discrete force operators, or an inaccurate curvature field. Balancing a constant curvature does not remove spatially varying curvature error. The curvature calculation and force balance are therefore separate requirements. The cell-integrated energy contribution uses the same momentum force: ```{math} :label: chap13_surface_tension_cell_work \int_{\Omega_P}\vec f_\sigma\cdot\vec u\,dV \simeq \Omega_P\vec f_{\sigma,P}\cdot\vec u_P ``` The force and velocity use consistent time levels. For a closed interface advected with the fluid, constant $\sigma$, and no contact-line work, the continuum relation to interface area $A_\Gamma$ is ```{math} :label: chap13_surface_tension_work \int_\Omega\vec f_\sigma\cdot\vec u\,dV =-\frac{d}{dt}\left(\sigma A_\Gamma\right) ``` Thus capillary work exchanges fluid and interfacial energy. Reusing the force does not alone guarantee exact discrete conservation of their sum. Phase change moves the interface relative to the fluid, so area change cannot then be identified with this work alone. ```{note} The optional {ref}`chap13_surface_diffusion` adds velocity diffusion along the interface. It remains in the converged momentum balance and damps physical motion as well as numerical motion. It does not correct curvature or, by itself, establish that a larger capillary timestep is accurate. ``` (chap13_wall_contact)= ### Wall contact The contact angle $\theta$ is measured through the color-one material where the interface meets a solid wall. With $\vec n_w$ pointing out of the fluid domain, a prescribed static angle imposes ```{math} :label: chap13_static_contact_angle \vec n\cdot\vec n_w=\cos\theta ``` Reversing the material pair replaces $\theta$ by $\pi-\theta$. Divergence methods impose the angle through the wall-normal flux; differential methods can also adjust wall-adjacent gradients. Geometric fits instead constrain the interface slope. Their constructions are given in {ref}`chap13_differential_wall` and {ref}`chap13_geometric_wall`. An uncorrected wall normal is not a prescribed 90-degree angle. A velocity-dependent contact-angle law is ```{math} :label: chap13_dynamic_contact_angle \theta=\operatorname{clip}\!\left( \sqrt[3]{\theta_{\mathrm{eq}}^3+b\,\mu U_{\mathrm{cl}}/\sigma}; \theta_{\min},\theta_{\max}\right) ``` Here $\theta_{\mathrm{eq}}$ is the prescribed equilibrium angle, and $\theta_{\min},\theta_{\max}$ bound the result. Angles are in radians, $b\ge0$ is a dimensionless coefficient, and $\mu$ is the adjacent-cell laminar viscosity. The signed speed $U_{\mathrm{cl}}$ is the adjacent-cell velocity relative to the wall, projected onto the wall-tangent direction normal to the contact line. Positive speed advances the color-one material. An unresolved tangent gives zero sampled speed. The law approximates contact-line speed by the adjacent fluid velocity and assumes a single equilibrium angle. A static angle prescribes the interface direction; contact-line motion also depends on the near-wall flow. (chap13_three_material_interfaces)= ### Three-material interfaces At a triple junction, a pair color also changes where the third material enters. Its gradient cannot represent only the selected pair's interface. The three-material model instead assigns a coefficient $\gamma_m$ to each material's entire boundary {ref}`[15] `. An A/B interface belongs to both boundaries, so $\gamma_A+\gamma_B=\sigma_{AB}$. Matching all three pair tensions gives ```{math} :label: chap13_three_tension_coefficients \gamma_A=\tfrac12(\sigma_{AB}+\sigma_{AC}-\sigma_{BC}) ``` The other coefficients follow by cyclic permutation of A, B, C. They are determined by the physical pair tensions, not additional material properties. Define the geometry colors by ```{math} :label: chap13_three_color_partition C_m=\frac{\alpha_m}{\alpha_A+\alpha_B+\alpha_C} \qquad \sum_m C_m=1 ``` Apply the same smoothing and nodal divergence to all three colors to obtain $\kappa_m$. Identical weights preserve the color sum and a material absent throughout the stencil. This normalization and smoothing do not change the transported material masses. The unsmoothed colors locate the force: ```{math} :label: chap13_three_material_force \vec f_\sigma=\sum_{m=A,B,C}\gamma_m\kappa_m\nabla C_m ``` If C is absent, $C_B=1-C_A$ and $\kappa_B=-\kappa_A$. The force reduces to ```{math} :label: chap13_three_binary_limit \vec f_\sigma =(\gamma_A+\gamma_B)\kappa_A\nabla C_A =\sigma_{AB}\kappa_A\nabla C_A ``` The surviving pair therefore retains its specified tension. Its discrete force, face pressure change, and energy work use the same construction as {ref}`chap13_capillary_balance`, summed over materials. This model requires positive pair tensions satisfying strict triangle inequalities, so all $\gamma_m>0$. Complete spreading and arbitrary four-or-more-material tensions require a different closure. At triple junctions, smoothing distributes the force over a finite region; the resulting angles depend on the curvature and force approximations {ref}`[16] `. Static wall angles must also satisfy a common wall-energy balance. Two independent pair angles determine the third; they cannot all be chosen independently. The equilibrium junction angles and coupled wall treatment are derived in {ref}`chap13_three_material_details` for static wall contact without phase change. (chap13_phase_change)= ## Interface Heat Transfer and Phase Change Phase change transfers mass between the liquid and vapor species of a substance. The transfer model supplies sources to the species and energy balances; the resulting composition and thermodynamic state determine the volume change seen by pressure correction. CMIM distinguishes transfer across a resolved interface from conversion within a cell: | Model | Rate determined by | Source location | | --- | --- | --- | | Kunkelmann | Heat supply at saturation temperature | Redistributed onto the two sides of an interface | | Schrage | Molecular kinetics and heat supply | Redistributed onto the two sides of an interface | | Zwart, Merkle, Schnerr--Sauer | Pressure departure from saturation | Equal and opposite species sources in one cell | The surface derivation below considers one pure liquid and its mapped vapor species. Schrage also permits a gas containing other species. Each transfer channel carries one constituent; contributions from several channels add in the governing equations. ### Interface heat balance Let $\Gamma$ denote the interface and $j''$ the mass transferred per area per time, in $\mathrm{kg\,m^{-2}\,s^{-1}}$, positive for evaporation and negative for condensation. The interface has temperature $T_\Gamma$ and no stored mass or energy. The specific enthalpy required for conversion is ```{math} :label: chap13_kunkelmann_phase_enthalpy L(p,T_\Gamma)=h_v(p,T_\Gamma)-h_l(p,T_\Gamma) ``` Here $h_l$ and $h_v$ are the mapped liquid and vapor constituent enthalpies, with compatible reference energies. For a gas mixture, $h_v$ is the vapor constituent's enthalpy, not the gas-mixture enthalpy. Both are evaluated at the local mechanical pressure $p$ and interface temperature. The surface models require $L>0$. Define $q_l''$ and $q_v''$ as the heat fluxes entering the interface from the liquid and gas sides. Neglecting kinetic-energy changes and viscous work in the interfacial heat balance gives the Stefan condition: ```{math} :label: chap13_kunkelmann_stefan_balance j''L=q_l''+q_v'' ``` Positive net heat supply causes evaporation; negative supply causes condensation. The two sides can also exchange heat with no mass transfer when $q_l''=-q_v''$. The interface is reconstructed at the half-value contour of the liquid/vapor pair color. Temperatures are sampled from cells sufficiently close to pure liquid or pure gas material on each side; the gas may contain carrier species. For side $m\in\{l,v\}$, write its heat supply as ```{math} :label: chap13_surface_thermal_coefficients q_m''=G_m(T_m-T_\Gamma) ``` $G_m$ is conductance per interface area, in $\mathrm{W\,m^{-2}\,K^{-1}}$, and $T_m$ is the effective sample temperature. One sample at normal distance $d_m>0$ gives $G_m=\lambda_m/d_m$. Multiple samples and propagation through a bounded band of mixed cells give the weighted construction in {ref}`chap13_interface_heat_sampling`. $T_\Gamma$ is a local interface closure, not an additional transported field. CMIM still solves one temperature per cell. The neighboring samples supply the two thermal gradients that cannot be recovered from a mixed cell's common temperature alone. The conductivity includes the material molecular value and the bulk eddy-diffusivity contribution when turbulence is present. This approximation does not model suppression of turbulence at the interface or an interfacial wall function. ### Surface transfer models The Kunkelmann model {ref}`[7] ` assumes local saturation equilibrium and negligible kinetic resistance: ```{math} :label: chap13_kunkelmann_interface_state T_\Gamma=T_{sat}(p) ``` Substituting the sampled heat supplies into the Stefan condition determines the rate: ```{math} :label: chap13_kunkelmann_sampled_stefan_flux j''=\frac{G_l(T_l-T_\Gamma)+G_v(T_v-T_\Gamma)}{L(p,T_\Gamma)} ``` This closure uses a pure liquid and pure vapor. Antoine or integrated Clausius--Clapeyron correlations supply the saturation relation; their assumptions and validity ranges are given in {ref}`chap13_saturation_relations`. The common-pressure approximation omits saturation shifts due to capillary pressure and vapor recoil. Schrage {ref}`[10] ` includes molecular transfer kinetics. At a trial interface temperature $T$, its signed mass flux is ```{math} :label: chap13_schrage_mass_flux j_S''(T)=\frac{2\chi}{2-\chi}\frac{1}{\sqrt{2\pi R_v}} \left[\frac{p_{sat}(T)}{\sqrt{T}}-\frac{p_v}{\sqrt{T_{inc}}}\right] ``` $\chi\in(0,1]$ is the accommodation probability and $R_v=R_u/M_v$ the vapor species gas constant. The terms represent molecules leaving the liquid and incident vapor molecules. The incident temperature $T_{inc}$ is approximated either by the sampled gas temperature $T_v$ or by the trial interface temperature $T$. For a gas mixture, the incident vapor pressure is ```{math} :label: chap13_schrage_vapor_partial_pressure p_v=p\langle X_v\rangle \qquad X_v=\frac{\widetilde Y_v/M_v}{\sum_{s\in\mathcal K_{gas}}\widetilde Y_s/M_s} ``` The mass fractions belong to the gas material. The brackets denote normalized face-area sampling on that side, without conductivity weighting. Pure vapor gives $p_v=p$. This is an ideal-gas partial-pressure approximation, even when the constituent enthalpies use a real-fluid EOS; fugacity corrections are excluded. Equating the kinetic rate to the heat-balance rate gives one scalar equation for $T_\Gamma$: ```{math} :label: chap13_schrage_temperature_residual G_l(T_l-T_\Gamma)+G_v(T_v-T_\Gamma) -j_S''(T_\Gamma)L(p,T_\Gamma)=0 ``` The enthalpies are reevaluated at each trial temperature and the common mechanical pressure. A bracketed solve requires a root within the shared validity range of the saturation relation and material EOS models. Both temperature resolution and the heat-balance residual must converge. The Stefan condition then supplies the mass flux used by the source construction. (chap13_surface_transfer_coupling)= ### Conservative coupling A patch of area $A_\Gamma$ transfers mass at rate $j''A_\Gamma$. The Hardt--Wondra procedure {ref}`[8] ` spreads this rate onto pure-material cells on the two sides. It determines where the source is applied, while Kunkelmann or Schrage determines how much mass transfers. Evaporation and condensation are redistributed separately to retain both when their net rate is zero. The smoothing length is numerical and controls the width of the source distribution. Let $S_l^{pc}$ and $S_v^{pc}$ be the mapped species sources per cell volume. With $\dot M=\int_\Gamma j''\,dA$, normalization gives ```{math} :label: chap13_phase_change_source_closure \int_\Omega S_l^{pc}\,dV=-\dot M \qquad \int_\Omega S_v^{pc}\,dV=+\dot M ``` Other species receive no source from this channel. Material sources follow by summing their species sources. Because removal and addition occur in different cells, $\dot m_{tot}^{pc}=\sum_m S_m^{pc}$ can be nonzero locally; its domain integral is zero. This local source remains in mixture continuity. The smoothing and normalization are derived in {ref}`chap13_mass_redistribution`. Energy transfer has three contributions: carried liquid enthalpy, carried vapor enthalpy, and heat exchanged with the thermal sample cells. Carried enthalpy uses the mass-redistribution weights. Each sampled heat flux is applied with opposite sign to its sample cell, using the same thermal sampling weights, as derived in {ref}`chap13_interface_heat_sampling`. Denote the sum of these energy sources, including carried kinetic energy, by $\widehat S_E^{pc}$. Its integral satisfies ```{math} :label: chap13_phase_change_energy_closure \int_\Omega\widehat S_E^{pc}\,dV =\int_\Gamma j''(h_v-h_l)\,dA -\int_\Gamma(q_l''+q_v'')\,dA =0 ``` The carried kinetic energies cancel between the two sides. The Stefan condition balances the enthalpy change with the sampled heat supply. The two heat contributions remain necessary even when $j''=0$. ```{important} The one-sided heat treatment replaces ordinary Fourier conduction on links crossing the reconstructed interface. Same-side conduction remains active. Latent enthalpy is already included in the transferred species enthalpies; no additional latent-heat source is added. ``` The mass sources also change the occupied material volume. At fixed thermodynamic state, their contribution is ```{math} :label: chap13_phase_change_material_volume_source S_V^{pc}=\sum_m\frac{S_m^{pc}}{\rho_m} ``` For constant liquid and vapor densities, its integral is $\dot M(1/\rho_v-1/\rho_l)$. Pressure correction accounts for this expansion through the predicted species inventories and volume-closure error, without an additional independent pressure source. The surface model sets $\vec S_u^{pc}=\vec 0$, omitting explicit momentum transfer between the separated source cells. Gravitational energy and continuity-subtracted energy assembly are treated in {ref}`chap13_phase_change_energy_assembly`. Several separated pure-liquid interfaces may feed different vapor species into a common gas through Schrage channels. Each channel retains its own geometry, thermal samples, and redistribution; their sources add. A crossing Fourier connection is removed once. This model excludes three-material cut cells and multicomponent liquid equilibrium. ### Volumetric cavitation Volumetric models represent unresolved liquid/vapor conversion driven by pressure departure from saturation. Each pair identifies one substance by liquid and vapor species, which may belong to the same material or to different materials. The rate laws give evaporation and condensation per volume occupied by the participating pair. Multiplying their difference by that pair's cell volume fraction $\phi$ gives the transfer per cell volume: ```{math} :label: chap13_cavitation_cell_volume_rate \dot m^{vol}=\phi\left(\widehat{\dot m}_{evap} -\widehat{\dot m}_{cond}\right) ``` The rates on the right are nonnegative; positive $\dot m^{vol}$ means evaporation. Zwart uses a prescribed bubble radius and nuclei fraction; Merkle uses pressure departure and a reference flow time; Schnerr--Sauer relates bubble radius to vapor content. The pair fractions and rate laws are given in {ref}`chap13_cavitation_rates`. The sources are collocated: ```{math} :label: chap13_volumetric_phase_change_sources S_l^{pc}=-\dot m^{vol} \qquad S_v^{pc}=+\dot m^{vol} ``` Mixture mass is therefore conserved in each cell. If the species belong to one material, its mass source also vanishes, although its composition and volume change. At fixed $p,T$, the pair's volume source is ```{math} :label: chap13_cavitation_volume_source S_V^{pc}=\dot m^{vol}\left(\frac{1}{\rho_v}-\frac{1}{\rho_l}\right) ``` Here $\rho_l$ and $\rho_v$ are pure-constituent EOS densities. The pressure/EOS closure accounts for this change through the updated species inventories. Conversion has zero conservative mixture-energy source; the enthalpy/EOS closure determines the temperature as composition changes. Compatible liquid/vapor enthalpy references remain essential. All pairs share cell pressure and temperature, including the influence of nonparticipating species. These models exclude heat-limited interface transfer, gas dissolution, and multicomponent phase equilibrium. Molecular and turbulent species diffusion are neglected in a material containing both species of a cavitation pair. ### Nonlinear treatment Surface geometry is reconstructed from the accepted beginning-of-step volume fractions and held fixed during the nonlinear iterations: ```{math} :label: chap13_phase_change_accepted_geometry \Gamma^{k}=\Gamma(\vec\alpha^{n}) ``` Pure-material support for sampling and redistribution is also fixed, avoiding changes of support from unconverged volume fractions. Pressure, temperature, conductivity, interface thermodynamics, rates, and redistribution weights update with the current flow iterate. Geometry is rebuilt for the next step; this first-order treatment requires the timestep to resolve interface motion. For volumetric cavitation, the pair volume and liquid/vapor fractions are held at their accepted values during the step. Pressure, temperature, saturation pressure, and pure-constituent densities update during iteration. The current pure-constituent density is needed even when that phase is absent. Each energy solve holds interface temperature and thermal conductance fixed and linearizes the heat withdrawal with respect to the sample-cell enthalpy. Its diagonal and matching source cancel at convergence; the derivation is in {ref}`chap13_phase_change_energy_assembly`. Mass-transfer rates and the interface-temperature solve are not included in this local derivative. The coupling is therefore segregated, with no fully coupled mass/energy Jacobian. A source-only step requires a timestep small enough to leave nonnegative species inventories. (chap13_outer_iteration)= (chap13_segregated_nonlinear_algorithm)= ## Segregated Nonlinear Algorithm Each timestep uses a segregated Picard iteration to advance species densities, velocity, and total enthalpy. The completed species predictor feeds the energy solve, giving this coupling a block Gauss--Seidel ordering. Transport coefficients and face reconstructions remain based on the incoming iterate $k$. The accepted state at time $n$ supplies the initial iterate and time histories. With split chemistry, reactions are integrated once before the loop to supply the reacted species history; see {ref}`chap13_chemistry_integration`. These histories remain fixed during iteration, as does the surface-transfer geometry described in {ref}`chap13_phase_change`. ```{figure} media/cmim_nonlinear_iteration.svg :alt: One nonlinear sweep. Species determines material masses, composition, and mixture density before the enthalpy solve. Momentum also uses that density when predictor gravity is enabled. Thermal closure precedes pressure correction and is repeated afterward. An unconverged state returns to the next sweep. :name: chap13_nonlinear_iteration_flowchart :width: 100% :align: center Dependencies within one nonlinear sweep. The dashed arrow is the optional use of predicted mixture density in the paired momentum and face-flux gravity terms. ``` ### Predictor updates Assemble species transport and its sources from the incoming iterate, using the linearization in {ref}`chap13_species_discretization`. Each linear solve holds its coefficients and source fixed. The completed species predictor determines material masses, mixture density, and composition: ```{math} :label: chap13_predictor_mass_sums q_m^{*,k}=\sum_{s\in\mathcal K_m}z_{ms}^{*,k} \qquad \rho^{*,k}=\sum_m q_m^{*,k} \qquad \widetilde Y_s^{(m),*,k}=\frac{z_{ms}^{*,k}}{q_m^{*,k}} \quad(q_m^{*,k}>0) ``` Mixture density is therefore available before energy or pressure correction, without an EOS evaluation. Intrinsic material densities $\rho_m$ require the thermal closures below. The energy solve accounts for the composition change at incoming pressure and temperature, together with an enthalpy-weighted species-residual correction. This couples the energy predictor to the new species state without changing the converged conservation equations. Its physical energy residual is retained for convergence checking. The correction uses incoming velocity; it does not depend on the momentum predictor. The derivation is in {ref}`chap13_material_energy_coupling`. The momentum predictor supplies velocity and the momentum-interpolated face volume flux. When gravity uses predicted density, both its cell force and paired face-flux contribution use $\rho^{*,k}$. This adds a species-to-momentum dependency. Other momentum density coefficients remain at the incoming iterate. Momentum and energy transport use the incoming mass flux in Eq. {eq}`chap13_incoming_picard_mass_flux`. ### Pressure and thermal closure Before pressure correction, recover temperature from Eq. {eq}`chap13_temperature_from_enthalpy` using the predicted species, total enthalpy, and velocity at pressure $p^k$. Evaluate the intrinsic material densities $\rho_m^{\mathrm{pred},k}$, the material-volume error, and its fixed-enthalpy pressure response. Solve the pressure-correction equation {eq}`chap13_fixed_q_pressure_row` using this volume error and the difference between the momentum and species predictor volume fluxes. Correct pressure, velocity, and face volume flux while holding $z_{ms}^{*,k}$ and $h_t^{*,k}$ fixed. At the corrected pressure, recover temperature again with static enthalpy $h_t^{*,k}-\tfrac12\lVert\vec u^{\mathrm c,k}\rVert^2$. Reevaluate the intrinsic material densities and derive $\alpha_m^{k+1}=q_m^{*,k}/\rho_m^{k+1}$. Mixture density remains $\rho^{*,k}=\sum_m q_m^{*,k}$. The corrected face flux affects species transport on the next nonlinear sweep. (chap13_pressure_transport_checks)= ### Nonlinear convergence Equation residuals, species and volume-fraction updates, material-volume closure, and corrected-flux continuity must satisfy their tolerances. The state must also satisfy the EOS, enthalpy closure, and volume-fraction bounds. The conservative species-update check is ```{math} :label: chap13_species_update \max_P \frac{\displaystyle\sum_m\sum_{s\in\mathcal K_m} \left|z_{ms,P}^{*,k}-z_{ms,P}^{k}\right|} {\max\!\left(\rho_P^{*,k},\rho_P^{k}\right)} \le\tau_z ``` Here $\tau_z$ is the species-update tolerance. Scaling by mixture density avoids division by a vanishing species or material mass. The continuity check uses the summed species flux at the corrected volume flow, with incoming donor and reconstruction coefficients held fixed: ```{math} :label: chap13_corrected_q_mass_flux \dot m_f^{\mathrm c,k} =\sum_m F_{q,m,f}^{*,k} +\left(\dot\Omega_f^{\mathrm c,k}-\dot\Omega_f^k\right) \sum_m\frac{\partial F_{q,m,f}}{\partial\dot\Omega_f} ``` A corrected flow that introduces an absent material also requires the appearance check in Eq. {eq}`chap13_pressure_appearance_update`. If a check is unsatisfied, repeat with the corrected state as iterate $k+1$. ```{include} face_reconstruction.inc ``` ```{include} species_diffusion.inc ``` (chap13_chemistry_integration)= ## Chemistry Integration For unsplit chemistry, the source in the species predictor is $\alpha_m\dot\omega_s^{(m)}=(q_m/\rho_m)\dot\omega_s^{(m)}$. Since $\widetilde Y_s^{(m)}=z_{ms}/q_m$, its destruction coefficient in Eq. {eq}`chap13_species_linearized_source` is ```{math} :label: chap13_species_chemistry_sink \ell_{ms}^k =\max\!\left(0,-\frac{1}{\rho_m^k} \frac{\partial\dot\omega_s^{(m)}}{\partial\widetilde Y_s^{(m)}}\bigg|_k\right) ``` The derivative holds $q_m$, $\rho_m$, pressure, temperature, and the other species fractions fixed. The coefficient has units of $\mathrm{s^{-1}}$; it is zero for absent materials and wholly explicit sources. Split chemistry instead integrates the cell reactions over $\Delta t$ from the accepted state, at fixed volume and internal energy. The transport step uses the reaction endpoint $z_{ms}^{\mathrm{chem}}$ as its BDF1 history: ```{math} :label: chap13_species_split_chemistry_balance \frac{\Omega_P}{\Delta t} \left(z_{ms,P}^{*,k}-z_{ms,P}^{\mathrm{chem}}\right) +\sum_{f\in\mathcal F_P}\sigma_{Pf} \left(F_{z,ms,f}^{*,k}+F_{ms,f}^{\mathrm{lin},*,k}\right) =\Omega_P S_{ms,P}^{\mathrm{tr,rem},k} ``` Here $S_{ms}^{\mathrm{tr,rem}}$ contains only transfer sources not included in the reaction integration. Chemistry is not added again as a source. Reactions within a material conserve its mass; reactions between materials may change their individual masses while preserving the total. The species reference energies contain reaction heat, so the total-energy history stays unchanged. This is first-order source splitting. ```{include} material_appearance_face_reference.inc ``` ```{include} transport_variants.inc ``` ```{include} curvature_methods.inc ``` ```{include} three_material_capillarity.inc ``` ```{include} surface_diffusion.inc ``` ```{include} phase_change_details.inc ``` ```{include} equilibria_and_balances.inc ``` ## References 1. (chap13_ref_thakur2022)= S. Thakur, J. Wright, and C. Neal, "An Algebraic Volume of Fluid Method Based on NVD and TVD Schemes Coupled to an All-Speed CFD Solver," AIAA SciTech Forum (2022), doi:10.2514/6.2022-1203. 2. (chap13_ref_brackbill1992)= J. U. Brackbill, D. B. Kothe, and C. Zemach, "A continuum method for modeling surface tension," J. Comput. Phys. 100, 335-354 (1992). 3. (chap13_ref_francois2006)= M. Francois et al., "A Balanced-Force Algorithm for Continuous and Sharp Interfacial Surface Tension Models Within a Volume Tracking Framework," J. Comput. Phys. 213, 141-173 (2006). 4. (chap13_ref_denner2014)= F. Denner and B. van Wachem, "Fully-Coupled Balanced-Force VOF Framework for Arbitrary Meshes with Least-Squares Curvature Evaluation from Volume Fractions," Numer. Heat Transfer B 65, 218-255 (2014). 5. (chap13_ref_demou2022)= A. D. Demou, N. Scapin, M. Pelanti, and L. Brandt, "A pressure-based diffuse interface method for low-Mach multiphase flows with mass transfer," J. Comput. Phys. 448, 110730 (2022), doi:10.1016/j.jcp.2021.110730. 6. (chap13_ref_clayton2026)= B. Clayton, J. McConnell, and Clell J. Solomon, "Analysis, thermodynamics, and a numeric solver for a pressure-temperature equilibrium closure of the four-equation model," arXiv:2606.27726 (2026). 7. (chap13_ref_kunkelmann2011)= C. Kunkelmann, *Numerical Modeling and Investigation of Boiling Phenomena*, doctoral dissertation, Technische Universität Darmstadt (2011), doi:10.26083/tuprints-00002731. 8. (chap13_ref_hardt2008)= S. Hardt and F. Wondra, "Evaporation model for interfacial flows based on a continuum-field representation of the source terms," J. Comput. Phys. 227, 5871--5895 (2008), doi:10.1016/j.jcp.2008.02.020. 9. (chap13_ref_zwart2004)= P. J. Zwart, A. G. Gerber, and T. Belamri, "A Two-Phase Flow Model for Predicting Cavitation Dynamics," Fifth International Conference on Multiphase Flow, Paper 152 (2004). 10. (chap13_ref_schrage1953)= R. W. Schrage, *A Theoretical Study of Interphase Mass Transfer*, Columbia University Press (1953). 11. (chap13_ref_jibben2019)= Z. Jibben, N. N. Carlson, and M. M. Francois, "A paraboloid fitting technique for calculating curvature from piecewise-linear interface reconstructions on 3D unstructured meshes," Computers & Mathematics with Applications 78, 643--653 (2019), doi:10.1016/j.camwa.2018.09.009. 12. (chap13_ref_bogner2016)= S. Bogner, U. Rüde, and J. Harting, "Curvature estimation from a volume of fluid indicator function for the simulation of surface tension and wetting with a free surface lattice Boltzmann method," Physical Review E 93, 043302 (2016), doi:10.1103/PhysRevE.93.043302. 13. (chap13_ref_merkle1998)= C. L. Merkle, J. Z. Feng, and P. E. O. Buelow, "Computational Modeling of the Dynamics of Sheet Cavitation," Third International Symposium on Cavitation, Grenoble (1998), 307--311. 14. (chap13_ref_schnerr2001)= G. H. Schnerr and J. Sauer, "Physical and Numerical Modeling of Unsteady Cavitation Dynamics," Fourth International Conference on Multiphase Flow (2001), [KIT record](https://publikationen.bibliothek.kit.edu/20722001). 15. (chap13_ref_dritselis2022)= C. Dritselis and G. Karapetsas, "Open-source finite volume solvers for multiphase (n-phase) flows involving either Newtonian or non-Newtonian complex fluids," arXiv:2203.09870 (2022), https://doi.org/10.48550/arXiv.2203.09870. 16. (chap13_ref_zhao2023)= C. Zhao, J. Maarek, S. M. Taleghani, and S. Zaleski, "A hybrid continuum surface force for the three-phase VOF method," arXiv:2310.12826 (2023), https://doi.org/10.48550/arXiv.2310.12826.