Mathematical Formulation & Numerics

This guide details the governing equations, spatial discretization, time-marching schemes, and projection steps implemented in the solver.

---

1. Scope & Implemented Physics

The solver is a collocated finite-volume method on unstructured hexahedral control volumes, supporting two distinct runtime paths: * Constant-Density Baseline: constant-density incompressible / low-Mach projection solver ($\nabla \cdot \mathbf{u} = 0$). * Active Variable-Density Path: variable-density low-Mach solver using Cantera/PelePhysics thermodynamics, conservative $\rho Y_k$ and $\rho h$ transport, and a conservative low-Mach divergence source.

---

2. Notation & Discretization Definitions

For a cell $c$: * $V_c$: Cell volume ($m^3$) * $f$: Face of cell $c$ * $A_f$: Face area ($m^2$) * $d_f$: Normal perpendicular distance used for face gradients ($m$) * $\mathbf{n}_f$: Face normal, owner-to-neighbor orientation * $\mathbf{n}_{f,c}$: Face normal reoriented outward from cell $c$ * $F_f$: Volumetric face flux, owner-to-neighbor orientation ($m^3/s$) * $F_{f,c}$: Volumetric face flux reoriented outward from cell $c$ * $\dot{m}_f$: Mass face flux, owner-to-neighbor orientation ($kg/s$) * $\dot{m}_{f,c}$: Mass face flux reoriented outward from cell $c$

The stored volumetric face flux is: $$F_f = \mathbf{u}_f \cdot \mathbf{n}_f A_f$$

The stored mass face flux is: $$\dot{m}_f = \rho_f F_f$$

For an upwind transported scalar $\psi$: $$\psi_f^{\text{up}} = \begin{cases} \psi_c, & F_{f,c} \ge 0 \\ \psi_{\text{other}}, & F_{f,c} < 0 \end{cases}$$

---

3. Density & Pressure Variables

The solver manages two distinct pressure and density scales: * Projection Pressure ($p$): Hydrodynamic pressure driving velocity projections ($\nabla \cdot \mathbf{u} = S_{\text{projection}}$). * Thermodynamic Pressure ($p_0$): Uniform background thermodynamic pressure (background_press) used for Cantera/PelePhysics state evaluations: $h(T, Y, p_0)$, $\rho(T, Y, p_0)$, etc. * Flow Density (transport%rho): Active flow density. Constant in baseline mode ($\rho = \rho_{\text{params}}$) and coupled to thermodynamics in variable-density mode ($\rho \leftarrow \rho_{\text{thermo}}$). * Thermodynamic Density ($\rho_{\text{thermo}}$): Density returned from state equations: $\rho_{\text{thermo}} = \rho(T, Y, p_0)$.

---

4. Governing Low-Mach Reacting Flow Equations

4.1 Momentum Conservation

$$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g}$$ where $\boldsymbol{\tau}$ is the viscous stress tensor: $$\boldsymbol{\tau} = \mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3} (\nabla \cdot \mathbf{u}) \mathbf{I} \right)$$ Kinematic viscosity is updated as $\nu = \mu / \rho$ when variable viscosity is active.

4.2 Species Transport

In variable-density mode: $$\frac{\partial (\rho Y_k)}{\partial t} + \nabla \cdot (\rho \mathbf{u} Y_k) = \nabla \cdot (\rho D_k \nabla Y_k) + \dot{\omega}_k \quad (k = 1, \dots, N_s)$$

4.3 Sensible Enthalpy Transport

$$\frac{\partial (\rho h)}{\partial t} + \nabla \cdot (\rho \mathbf{u} h) = \nabla \cdot (\lambda \nabla T) + \nabla \cdot \left( \rho \sum_{k=1}^{N_s} h_k D_k \nabla Y_k \right) + \dot{q}_{\text{chem}} - \nabla \cdot \mathbf{q}_{\text{rad}}$$ Sensible enthalpy is defined relative to reference temperature $T_{\text{ref}}$ (typically $298.15\text{ K}$) to prevent artificial heat release during non-reacting mixing: $$h_{\text{sens}}(T, Y, p_0) = h_{\text{abs}}(T, Y, p_0) - h_{\text{abs}}(T_{\text{ref}}, Y, p_0)$$

---

5. Fractional-Step Pressure Projection Method

5.1 Momentum Predictor Step

The predictor computes a tentative cell velocity $\mathbf{u}^*$ using an explicit Adams-Bashforth 2nd-order (AB2) time scheme: $$\mathbf{u}^* = \mathbf{u}^n + \Delta t \left( \frac{3}{2} \mathbf{R}_{\mathbf{u}}^n - \frac{1}{2} \mathbf{R}_{\mathbf{u}}^{n-1} \right)$$ where the momentum residual is: $$\mathbf{R}_{\mathbf{u},c} = -\frac{1}{V_c} \sum_{f \in \partial c} F_{f,c} \mathbf{u}_f^{\text{adv}} + \frac{1}{V_c} \sum_{f \in \partial c} \nu_f A_f \frac{\mathbf{u}_{\text{other}} - \mathbf{u}_c}{d_f} - \frac{1}{\rho_c} \nabla p_c + \mathbf{g}$$ The pressure gradient is calculated using the Gauss finite-volume reconstruction: $$\nabla p_c = \frac{1}{V_c} \sum_{f \in \partial c} p_f \mathbf{n}_{f,c} A_f$$

5.2 Poisson Solver for Pressure Correction

The corrected face fluxes satisfy the low-Mach divergence constraint: $$\sum_{f \in \partial c} F_{f,c}^{n+1} = S_c V_c$$ The face correction is: $$F_f^{n+1} = F_f^* - \Delta t \frac{A_f}{\rho_f d_f} \left( \phi_{\text{nb}} - \phi_{\text{owner}} \right)$$ The elliptic system $A\phi = b$ is formulated as: $$(A\phi)_c = \sum_{f \in \partial c} \frac{A_f}{\rho_f d_f} \left( \phi_c - \phi_{\text{other}} \right)$$ $$b_c = -\frac{1}{\Delta t} \left[ \sum_{f \in \partial c} F_{f,c}^* - S_c V_c \right]$$

---

6. Preconditioned Conjugate Gradient (PCG) Solver

The Poisson system is solved iteratively using a parallel PCG method implemented in src/numerics/mod_linear_solver.f90:

1. Diagonal Jacobi Preconditioner ($M$): $$M_{c,c} = A_{c,c}^{-1} = \left[ \sum_{f \in \partial c} \frac{A_f}{\rho_f d_f} \right]^{-1}$$
2. Iteration Loop: For an initial guess $\phi^0$:
3. Compute initial residual: $\mathbf{r}^0 = \mathbf{b} - A\boldsymbol{\phi}^0$, apply preconditioner: $\mathbf{z}^0 = M\mathbf{r}^0$, set search direction: $\mathbf{p}^0 = \mathbf{z}^0$.
4. For each iteration $k$:
   • Compute matrix-vector product: $\mathbf{w}^k = A\mathbf{p}^k$
   • Compute step length: $\alpha^k = (\mathbf{r}^k, \mathbf{z}^k) / (\mathbf{p}^k, \mathbf{w}^k)$
   • Update solution: $\boldsymbol{\phi}^{k+1} = \boldsymbol{\phi}^k + \alpha^k \mathbf{p}^k$
   • Update residual: $\mathbf{r}^{k+1} = \mathbf{r}^k - \alpha^k \mathbf{w}^k$
   • Check convergence: if $|\mathbf{r}^{k+1}\|_2 \le \epsilon_{\text{atol}} + \epsilon_{\text{rtol}} |\mathbf{b}\|_2$, exit.
   • Apply preconditioner: $\mathbf{z}^{k+1} = M\mathbf{r}^{k+1}$
   • Compute update factor: $\beta^k = (\mathbf{r}^{k+1}, \mathbf{z}^{k+1}) / (\mathbf{r}^k, \mathbf{z}^k)$
   • Update search direction: $\mathbf{p}^{k+1} = \mathbf{z}^{k+1} + \beta^k \mathbf{p}^k$
5. Nullspace Projection (Pressure Pin): For pure-Neumann boundaries, the constant pressure nullspace is removed by projecting a zero-mean gauge: $$\phi_c \leftarrow \phi_c - \frac{1}{\sum V_i} \sum V_i \phi_i$$

---

7. High-Order Spatial Reconstruction Schemes

To reduce numerical diffusion while enforcing physical boundedness, the solver supports high-order convection schemes:

7.1 Quadratic Upstream Interpolation (QUICK)

For a face $f$ with flow from cell $C$ to cell $D$: $$\psi_f = \frac{1}{2}(\psi_C + \psi_D) - \frac{1}{8}(\psi_U - 2\psi_C + \psi_D)$$ The virtual upstream node $U$ is reconstructed dynamically on unstructured grids: $$\psi_U = 2\psi_C - \psi_D + 2 \nabla\psi_C \cdot (\mathbf{r}_C - \mathbf{r}_D)$$

7.2 Hybrid Linear/Parabolic Approximation (HLPA) Limiter

To prevent unphysical oscillations, QUICK is combined with a low-order upwind fallback using the Normalized Variable Diagram (NVD): $$\tilde{\psi}_C = \frac{\psi_C - \psi_U}{\psi_D - \psi_U}$$ $$\tilde{\psi}_f = \begin{cases} \tilde{\psi}_C + \frac{1}{2}(1 - \tilde{\psi}_C) & \text{if } 0 < \tilde{\psi}_C < 1 \\ \tilde{\psi}_C & \text{otherwise} \end{cases}$$ $$\psi_f = \psi_U + \tilde{\psi}_f (\psi_D - \psi_U)$$

---

8. Stiff Chemistry Integration & Adaptive Subcycling

Decoupled from transport via a fractional-step operator-splitting approach, reactions are integrated using CVODE.

8.1 Adaptive Subcycling

The global flow step $\Delta t$ is divided into a cell-local number of sub-slices $N_{\text{sub}}$: $$\Delta t_{\text{sub}} = \frac{\Delta t}{N_{\text{sub}}}$$ $$N_{\text{sub}} = \max\left(1, \min\left(N_{\text{max}}, \left\lceil \frac{|T^n - T^{n-1}|}{\Delta T_{\text{threshold}}} \right\rceil\right)\right)$$ where $N_{\text{max}}$ is max_chemistry_subcycles and $\Delta T_{\text{threshold}}$ is subcycling_dT_threshold ($10\text{ K}$).

8.2 Dominant-Species Simplex Projection

To preserve mass conservation without triggering numerical stiffness in trace intermediate radicals, mass fractions are normalized by modifying only the dominant majority bath species (typically $N_2$): $$\delta = 1.0 - \sum_{k=1}^{N_s} Y_k \implies Y_{\text{bath}} \leftarrow Y_{\text{bath}} + \delta$$

8.3 Continuous $\mathcal{C}^1$ Cubic Spline Temperature Clipping

To avoid non-differentiable step transitions that break Newton-solver Jacobian evaluations, temperature clipping is smoothed using a cubic spline over window $[T_{\text{clip}}, T_{\text{clip}} + \Delta T_{\text{clip}}]$: $$T_{\text{clipped}} = \begin{cases} T & T \ge T_{\text{clip}} + \Delta T_{\text{clip}} \\ T_{\text{clip}} + a(T - T_{\text{clip}})^2 + b(T - T_{\text{clip}})^3 & T_{\text{clip}} < T < T_{\text{clip}} + \Delta T_{\text{clip}} \\ T_{\text{clip}} & T \le T_{\text{clip}} \end{cases}$$ $$a = \frac{3}{\Delta T_{\text{clip}}}, \quad b = -\frac{2}{(\Delta T_{\text{clip}})^2}$$

---

9. Dynamic Timestep Control, Rejection & Damping

9.1 Timestep Rejection & Rollback

If a physical check fails, the timestep is rejected, the state is rolled back to time level $n$, and retried with $\Delta t_{\text{new}} = \Delta t_{\text{old}} \times \text{step\_reject\_cut\_factor}$. Rejections are triggered if: * NaN/Inf states occur. * Pressure Poisson CG iterations exceed pressure_max_iter. * Local CFL exceeds stability limits. * Density ratios fall below safety thresholds.

9.2 Chemical Source Damping

To prevent pressure shocks in low-Mach projection source terms due to massive local energy release, chemical source terms can be relaxed: $$S_{\text{chem}}^{\text{effective}} = S_{\text{chem}}^{\text{old}} + \alpha_{\text{relax}} \left( S_{\text{chem}}^{\text{new}} - S_{\text{chem}}^{\text{old}} \right)$$ where $\alpha_{\text{relax}}$ is chemistry_source_relaxation.