8 Discretization of the Dispersed-Phase Equations

8.1 Finite-volume form

For a dispersed conserved quantity vector \(\mathbf{U}_d\), the cell balance is \[\frac{d}{dt}\int_{\Omega_P}\mathbf{U}_d\,d\Omega +\sum_{f\in\partial P}\mathbf{F}_{d,f}A_f =\Omega_P\mathbf{S}_{d,P}. \](8.1) The active dispersed block contains mass, momentum, a thermal variable, and optionally interfacial area. Source terms include interphase momentum/heat transfer, body forces, and IATE breakup/coalescence contributions.

8.2 Primitive-to-transported-variable dependence

The finite-volume residual is assembled from primitive variables, while the stored/transported particle quantities are nonlinear functions of them. In particular, \[\alpha=\frac{\rho_d}{\rho_{p,m}(T_d)}, \qquad d_d=\frac{6\rho_d}{\rho_{p,m}(T_d)a_i}.\] When IATE is active, a perturbation of \(a_i\) therefore changes the diameter, and a perturbation of \(T_d\) can change the diameter through the condensed-material density. AD retains these dependencies in the source and flux Jacobians.

8.3 Regularized normal subsystem

With packing friction pressure active, the one-dimensional mass/momentum subsystem normal to a face is \[\mathbf{U}_{d,n}=\begin{bmatrix}\rho_d\\\rho_du_{n,d}\end{bmatrix}, \qquad \mathbf{F}_{d,n}=\begin{bmatrix} \rho_du_{n,d}\\ \rho_du_{n,d}^{2}+p_{fr} \end{bmatrix}.\] The signal estimate is \[s_d=|u_{n,d}|+c_d, \](8.2) where \(c_d\) is given by Eq. (7.18). In the dilute pressureless limit \(c_d\rightarrow0\).

8.4 Rusanov-type particle flux

For a left/right state, the local Lax–Friedrichs/Rusanov particle flux is \[\mathbf{F}_{d,f}^{*} =\frac12\left(\mathbf{F}_{d,L}+\mathbf{F}_{d,R}\right) -\frac12s_{\max}\left(\mathbf{U}_{d,R}-\mathbf{U}_{d,L}\right), \](8.3) with \[s_{\max}=\max\left( |u_{n,d,L}|+c_{d,L}, |u_{n,d,R}|+c_{d,R} \right).\] The regularization therefore adds numerical diffusion only where the friction-pressure model produces a nonzero \(c_d\).

8.5 Interfacial-area convective flux

When IATE is enabled, interfacial area is advected with the dispersed transport velocity. A first-order upwind representation is \[F_{a_i,f}=a_{i,\mathrm{up}}u_{n,d,f}. \](8.4) The IATE row receives the local volumetric source \[S_{IATE}=S_{RC}+S_{RT}+S_{TI}+S_G.\] Because \(d_d\propto 1/a_i\), positivity and boundedness of the transported area are physically significant, not merely scalar-transport niceties.

8.6 Receiver-side packing limiter

A conservative face flux can create an inadmissible receiving-cell state if the transported dispersed mass exceeds local packing. For a trial transfer \(\Delta m_f\), CMPS limits the accepted transfer as \[\Delta m_f^{\mathrm{acc}} =\theta_f\Delta m_f, \qquad 0\le\theta_f\le1, \](8.5) where \(\theta_f\) is selected so that \[\rho_d^{n+1}\le\rho_{d,\max} =\alpha_{\max}\rho_{p,m}(T_d). \](8.6) The accepted factor must be applied consistently to the dispersed quantities transported by that face so that mass, momentum, thermal transport, and area transport do not correspond to different effective particle fluxes.

8.7 Cell-volume scaling of IATE sources

The continuum source \(S_{IATE}\) has units \(\mathrm{m}^{-1}\mathrm{s}^{-1}\). Its finite-volume contribution to the cell residual is therefore \[R_{a_i,P}^{src}=-\Omega_P S_{IATE,P} \](8.7) for the residual convention in which fluxes minus sources define \(R\). The global Newton right-hand side remains \(B=-R\) according to Chapter 3.

This distinction prevents a common sign error. A physical breakup source is positive in the differential IATE, but with a residual written as “transport minus source”, its direct residual contribution is negative. The assembled Newton correction nevertheless follows the common CMPS \(A\,\delta q=B\) convention.

8.8 Source Jacobians and local cross-coupling

For a generic local source \(\mathbf{S}_{dg}(\mathbf{q}_g,\mathbf{q}_d)\), AD produces \[\frac{\partial\mathbf{S}_{dg}}{\partial\mathbf{q}_g}, \qquad \frac{\partial\mathbf{S}_{dg}}{\partial\mathbf{q}_d}.\] For the IATE row specifically, \[\frac{\partial R_{a_i,P}^{src}}{\partial q_j} =-\Omega_P\frac{\partial S_{IATE,P}}{\partial q_j}. \](8.8) The derivatives can include the chain \[q_j\rightarrow \{\alpha,d_d,\mathrm{Re}_d,\epsilon_t,\lambda_B,\ldots\} \rightarrow S_{IATE},\] so a fully implicit IATE is substantially more coupled than an advected passive scalar.

For example, the diameter derivatives for constant \(\rho_{p,m}\) are \[\frac{\partial d_d}{\partial\rho_d}=\frac{d_d}{\rho_d}, \qquad \frac{\partial d_d}{\partial a_i}=-\frac{d_d}{a_i}. \](8.9) If \(\rho_{p,m}=\rho_{p,m}(T_d)\), an additional derivative appears, \[\frac{\partial d_d}{\partial T_d} =-\frac{d_d}{\rho_{p,m}} \frac{d\rho_{p,m}}{dT_d}. \](8.10) These relations illustrate why source linearization should be obtained from the same active expressions used for the residual rather than from manually simplified partial derivatives.

8.9 Turbulent-impact source discretization

For the turbulent-impact model of Eq. (7.27), the cell source is evaluated from the local active state, including \[\epsilon_t=0.09k\omega, \qquad d_d=\frac{6\rho_d}{\rho_{p,m}a_i}, \qquad \lambda_B =\exp\!\left[-\frac{K_B\sigma} {\rho_gd_d^{5/3}\epsilon_t^{2/3}}\right].\] The implementation protects small denominators and the remaining-packing expression. This is numerically necessary because an apparently benign source can become extremely stiff when \(a_i\rightarrow0\), \(d_d\rightarrow0\), turbulence vanishes, or the packing margin collapses.

The same protected algebra must be differentiated by AD. Applying a protection only after derivative evaluation would make the residual and Jacobian inconsistent.

8.10 Source time scales and stiffness

A local IATE source time scale can be estimated as \[\tau_{IATE}=\frac{a_i}{|S_{IATE}|+S_{floor}}, \](8.11) where \(S_{floor}\) is a numerical floor used only for estimating the scale. If \(\tau_{IATE}\) becomes much smaller than the carrier or particle convective time, explicit source treatment would require severe time-step restriction. This is one reason the CMPS formulation includes source derivatives in the coupled implicit system.

Analogously, the drag and heat-transfer response times \(\tau_v\) and \(\tau_T\) from Chapter 7 indicate when carrier/dispersed momentum and temperature equations become strongly coupled. The nonlinear solve is most demanding when several of these time scales become simultaneously short.

8.11 Boundary treatment of dispersed variables

At an inlet where dispersed material is prescribed, the boundary state must specify enough information to reconstruct the active particle variables. In monodisperse mode this includes the prescribed representative diameter. In IATE mode a prescribed diameter \(d_{d,b}\) and dispersed mass concentration \(\rho_{d,b}\) imply \[a_{i,b}=\frac{6\rho_{d,b}} {\rho_{p,m}(T_{d,b})d_{d,b}}. \](8.12) Using an inconsistent independent \(a_i\) and \(d_d\) at the same boundary would violate the one-group closure.

Wall/outlet behavior depends on the selected dispersed boundary model. Any boundary construction used during implicit assembly must preserve positivity and must provide a differentiable state wherever its dependence is included in the Jacobian.

8.12 Conservation and model qualification

Breakup and coalescence redistribute interfacial area but do not create or destroy dispersed mass. Therefore the IATE source must not be inserted into the dispersed mass row. Momentum exchange should be equal and opposite between carrier and dispersed phases. Heat exchange should likewise be paired with opposite signs, subject to the definition of the transported thermal variables.

The August 2026 implementation audit found that the active drag-energy work partition does not yet provide the desired combined carrier-plus-particle energy balance for nonzero slip. This issue is independent of the IATE geometric closure and is recorded in Chapter 20 rather than hidden inside the dispersed discretization.

8.13 Current dispersed-phase convective flux family

The current dispersed-phase options expose four face-flux families: the standard donor/receiver flux, Rusanov, AUSM and HLLC. All four transport dispersed mass, momentum, total energy and, when IATE is active, interfacial area. Packing and minimum-density gates prevent a donor below the active dispersed-density floor from emitting particles and prevent transport into a receiver already at its maximum material packing density.

8.13.1 Standard donor/receiver flux

The standard path uses the signs of the left and right particle normal velocities. Same-direction motion selects the appropriate donor. Opposing motion toward the face permits contributions from both sides; motion away from the face gives zero convective particle flux. For a left-to-right donor, for example,

\[\dot m_d=\rho_{d,L}u_{n,L},\quad \mathbf F_m=\dot m_d\mathbf V_{d,L},\quad F_E=\dot m_dE_{d,L},\quad F_{a_i}=u_{n,L}a_{i,L}.\]

8.13.2 Rusanov flux

With frictional-pressure wave speeds \(c_{d,L},c_{d,R}\), the dissipation speed is

\[a_{max}=\max(|u_{n,L}|+c_{d,L},\ |u_{n,R}|+c_{d,R}).\]

The mass flux is

\[F_\rho=\frac12(\rho_Lu_{n,L}+\rho_Ru_{n,R})-\frac12a_{max}(\rho_R-\rho_L),\]

with analogous central-minus-jump terms for momentum, energy and interfacial area. The face friction pressure is the arithmetic mean and is added only to the momentum flux.

8.13.3 Dispersed AUSM flux

The dispersed AUSM path applies Mach/pressure splitting to the regularized particle system. Friction pressure supplies the pressure part of momentum flux, while energy is convected with the particle mass transport. The final mass-flux direction is passed through the same donor/receiver packing gate used by the other particle fluxes.

8.13.4 Dispersed HLLC flux

The HLLC path treats the artificial friction pressure as the pressure variable of the regularized particle hyperbolic system. Its bounding waves are

\[S_L=\min(u_{n,L}-c_{d,L},u_{n,R}-c_{d,R}),\qquad S_R=\max(u_{n,L}+c_{d,L},u_{n,R}+c_{d,R}),\]

and the contact wave is

\[S_M=\frac{p_{fr,R}-p_{fr,L}+\rho_Lu_{n,L}(S_L-u_{n,L})-\rho_Ru_{n,R}(S_R-u_{n,R})}{\rho_L(S_L-u_{n,L})-\rho_R(S_R-u_{n,R})}.\]

Left/right star densities, velocities, energy and interfacial area are reconstructed from \(S_M\); the final flux follows the usual four HLLC wave regions and is finally checked against the receiver packing limit.