CMPS Theory & User Reference Manual › Numerical Approach

5 Numerical Approach

5.1 Cell-centred finite-volume discretization

For each control volume \(P\), the spatial residual is assembled from numerical convective fluxes, viscous/diffusive fluxes, and cell sources, \[\mathbf{R}_P(\mathbf{q}) =\sum_{f\in\partial P} \left(\mathbf{F}^{c}_f-\mathbf{F}^{v}_f\right)A_f -\Omega_P\mathbf{S}_P. \](5.1) The unknown vector is cell-centred. Interior faces couple two cells and therefore generate both diagonal and neighbour Jacobian blocks. Boundary faces generate owner-cell blocks whose exact form depends on the boundary condition.

The steady nonlinear problem is \[\mathbf{R}(\mathbf{q})=\mathbf{0}.\] Linearization about iteration \(m\) gives \[\frac{\partial\mathbf{R}}{\partial\mathbf{q}}\bigg|_m \delta\mathbf{q} =-\mathbf{R}(\mathbf{q}^m).\] With the CMPS storage convention \(\mathbf{B}=-\mathbf{R}\) and \(\mathbf{A}=\partial\mathbf{R}/\partial\mathbf{q}\), this becomes Eq. (3.3).

5.2 Implicit pseudo-time stabilization

The steady solver may add a local pseudo-time mass matrix to improve nonlinear robustness, \[\left[ \frac{\Omega_P}{\Delta\tau_P}\boldsymbol{\Gamma}_P +\frac{\partial\mathbf{R}_P}{\partial\mathbf{q}} \right]\delta\mathbf{q} =-\mathbf{R}_P, \](5.2) where \[\boldsymbol{\Gamma}=\frac{\partial\mathbf{W}_{\tau}}{\partial\mathbf{q}}\] is the pseudo-time transformation between primitive corrections and the stored pseudo-state.

For compressible calculations, \(\boldsymbol{\Gamma}\) can be modified by local time-derivative preconditioning to reduce the disparity between convective and acoustic scales at low Mach numbers (Choi and Merkle 1993; Weiss and Smith 1995). For the constant-density incompressible regime, the pressure-row coefficient is \(1/c_{ac}^{2}\) as described by Eq. (4.2).

The pseudo time step is local and is selected from a characteristic spectral scale and the requested CFL number. Increasing CFL reduces the diagonal pseudo-time stabilization and moves the iteration toward a Newton solve; decreasing CFL increases diagonal dominance and generally improves robustness at the cost of more nonlinear iterations. Automatic CFL logic should therefore be interpreted as nonlinear continuation rather than as physical time integration.

5.2.1 Isentropic-Mach preconditioning accelerator

The low-Mach preconditioner uses a local reference velocity \(u_{\mathrm{ref}}\). If this scale becomes too small in a strong expansion, the preconditioned acoustic system can become more aggressive than the local thermodynamic state warrants. CMPS provides an isentropic-Mach accelerator that raises the reference-velocity floor from a total-to-static pressure estimate.

A reference total pressure \(p_{t,\infty}\) is obtained from the first applicable inlet-type boundary. For a far-field boundary, the source evaluates

\[p_{t,\infty}=\left(p_{\infty}+p_{\infty}^{\mathrm{EOS}}\right)\left(1+\frac{\gamma_{\infty}-1}{2}M_{\infty}^2\right)^{\gamma_{\infty}/(\gamma_{\infty}-1)}-p_{\infty}^{\mathrm{EOS}},\]

where \(p_{\infty}^{\mathrm{EOS}}=p_{\infty,\mathrm{stiff}}\) for a stiffened-gas material and zero otherwise. Mass-flow-inlet and solid-propellant-surface boundaries provide their stored reference total pressure directly, while a stagnation inlet provides its specified total pressure. If no suitable inlet is found, the reference value is zero and the isentropic limiter becomes inactive after pressure-ratio clipping.

Before the isentropic limiter is applied, the local reference scale is assembled from the acoustic floor, the largest local or neighboring velocity, pressure jumps, diffusion and, when selected, a global reference Mach number. Its structure can be summarized as

\[u_{\mathrm{ref},0}=\max\!\left(M_{r,\min}a,\ u_{\max},\ 2\sqrt{\frac{\Delta p_{\max}}{\rho_f}},\ u_d,\ M_{\mathrm{global}}a\right),\qquad u_d=\frac{\lambda_v}{L_c},\]

with the global term present only when global preconditioning is enabled. The isentropic pressure ratio is then

\[\Pi_{\mathrm{is}}=\max\!\left(\frac{p_{t,\infty}+p_{\mathrm{shift}}}{p+p_{\mathrm{shift}}},\ 1\right),\]

and the corresponding isentropic Mach estimate is

\[\boxed{M_{\mathrm{is}}=\sqrt{\frac{2}{\gamma-1}\left(\Pi_{\mathrm{is}}^{(\gamma-1)/\gamma}-1\right)}}.\]

The preconditioning reference velocity is raised according to

\[u_{\mathrm{ref}}\leftarrow\max\!\left(u_{\mathrm{ref},0},\ M_{\mathrm{is}}a\right).\]

For transient calculations the physical-time scale supplies an additional lower bound,

\[u_{\mathrm{ref}}\leftarrow\max\!\left(u_{\mathrm{ref}},\ \frac{L_c}{\pi\Delta t}\right),\]

and, when a finite acoustic speed exists, CMPS finally limits the reference velocity by

\[u_{\mathrm{ref}}\leftarrow\min\!\left(u_{\mathrm{ref}},a\right).\]

Stiffened-gas treatment. In the local isentropic limiter, \(p_{\mathrm{shift}}=p_{\infty,\mathrm{stiff}}\) only for a non-VOF stiffened-gas carrier. Otherwise \(p_{\mathrm{shift}}=0\). Hence the non-VOF stiffened-gas form is

\[\Pi_{\mathrm{is}}=\max\!\left(\frac{p_{t,\infty}+p_{\infty,\mathrm{stiff}}}{p+p_{\infty,\mathrm{stiff}}},1\right),\qquad M_{\mathrm{is}}=\sqrt{\frac{2}{\gamma-1}\left(\Pi_{\mathrm{is}}^{(\gamma-1)/\gamma}-1\right)}.\]

VOF simulation note. When homogeneous VOF is enabled, the current local limiter sets \(p_{\mathrm{shift}}=0\) even when a constituent uses a stiffened-gas EOS. The mixture acoustic speed remains the speed used to convert \(M_{\mathrm{is}}\) into the velocity floor.

5.3 5.3 Consistent residual linearization

CMPS forms the implicit face Jacobians from the same numerical expressions used in the finite-volume residual. For an interior face, the linearized flux is

5.4 Reconstruction and order of accuracy

A face state may be written generically as \[\phi_f^{L} =\phi_L+\Psi_L\nabla\phi_L\cdot(\mathbf{x}_f-\mathbf{x}_L),\] with an analogous expression for the right cell. \(\Psi\) denotes the active limiter/reconstruction control. First-order reconstruction sets the correction to zero and is the most dissipative, while higher-order reconstruction improves spatial accuracy but increases sensitivity to mesh quality and nonlinear overshoots.

For an implicit iteration, the reconstructed residual and the Jacobian must represent the same frozen reconstruction state used in that assembly. Deferred-correction or frozen-gradient strategies are valid only when their lagging is deliberate and consistently reflected in the nonlinear iteration.

5.5 Convective and pressure fluxes

The compressible path uses the selected Riemann/AUSM-family carrier flux. For AUSM⁺-up, a face flux may be represented schematically as \(\mathbf F_f^c=\dot m_f\boldsymbol\Phi_{up}+p_f\boldsymbol\Pi_n\). The pressure-difference contribution in mass flux and the velocity-difference contribution in pressure are included in both the residual and its implicit linearization, which strengthens pressure-velocity coupling at low Mach number.

The constant-density incompressible path uses the same general split structure with \(c_{ac}\) as the characteristic speed; its exact row fluxes are given by Eq. (4.4).

5.6 Viscous and diffusive terms

Viscous stresses, heat conduction, turbulence diffusion and species diffusion are assembled as face fluxes. A typical scalar contribution is \(F_{\phi,f}^{v}=\Gamma_{\phi,f}(\nabla\phi)_f\cdot\mathbf n_f\). On non-orthogonal meshes the face gradient includes the selected reconstruction and geometric correction. The same finite-volume sign convention is used in the implicit matrix, so stronger diffusion generally improves damping but also increases coupling and linear-system stiffness.

5.7 Species admissibility and coupled relaxation

The \(N-1\) species formulation requires \[Y_s\ge0, \qquad \sum_{s=1}^{N-1}Y_s\le1.\] The current coupled update applies a common per-cell limiter when an unconstrained increment would violate the independent-species bounds. Scaling the full cell correction rather than clipping each variable independently better preserves the direction of the coupled Newton correction. Additional physical limits, such as temperature bounds, are applied by the property/update path where configured.

5.8 Physical transient discretization

For transient calculations, the physical-time residual is Eq. (3.4). BDF1 is first-order accurate and requires one previous accepted physical state. BDF2 is second-order accurate and requires two historical states. The spatial residual at the new physical time is solved implicitly through pseudo-time iterations.

The dual-time nonlinear problem at physical step \(n+1\) is \[\mathbf{R}^{*}(\mathbf{q}^{n+1}) =\mathbf{R}^{\mathrm{space}}(\mathbf{q}^{n+1}) +\mathbf{R}^{\mathrm{time}}(\mathbf{q}^{n+1};\mathbf{W}^n,\mathbf{W}^{n-1}) =\mathbf{0}. \](5.4) An inner pseudo-time iteration solves \[\left[ \frac{\Omega_P}{\Delta\tau_P}\boldsymbol{\Gamma}_{\tau,P} +\frac{\partial\mathbf{R}^{*}}{\partial\mathbf{q}} \right]\delta\mathbf{q} =-\mathbf{R}^{*}. \](5.5) Only the accepted physical-time solution advances the BDF history.

For constant-density artificial-compressibility flow, the pressure row has no physical \(dp/dt\) term. The pressure pseudo-time term remains in \(\boldsymbol{\Gamma}_{\tau}\) to converge the divergence constraint within each physical step.

5.9 Residual monitoring

A raw residual component has physical units and depends on equation scaling, mesh area/volume, and state magnitude. CMPS therefore monitors scaled residual quantities in addition to the nonlinear update history. A generic normalized component may be written as \[\widehat R_j =\frac{\lVert R_j\rVert}{R_{j,\mathrm{ref}}},\] where the reference must be defined consistently over the run. Residual reduction alone is not a guarantee of physical admissibility; density, temperature, species, turbulence variables, and other model-specific bounds must also remain valid.

5.10 Relation between nonlinear and linear convergence

The Krylov solve in Eq. (3.3) is an inner problem. Driving it far below the accuracy justified by the current nonlinear linearization wastes work, whereas stopping it too early can destroy the intended Newton correction. The appropriate linear tolerance therefore depends on nonlinear stage, pseudo-time stabilization, and preconditioner quality. CMPS exposes the matrix-solver options separately from the nonlinear CFL/relaxation controls for this reason.

5.11 Coupled correction admissibility

The linear correction is not committed blindly. CMPS scales the complete cell correction when temperature, species or VOF bounds would be crossed, preserving the coupled direction rather than clipping individual solved variables independently. The exact temperature, species and VOF fraction-to-boundary equations are given in Nonlinear Correction, Relaxation and Admissibility.