4 Artificial Compressibility for Constant-Density Flow

The incompressible carrier solver is a dedicated constant-density formulation that reuses the coupled finite-volume infrastructure without treating pressure as an equation-of-state variable. This chapter documents the state, pseudo-time equation, inviscid flux, and the main implementation consequences.

4.1 Primitive and stored states

For dimension \(d\) and \(N\) species, the active primitive block has the form \[\mathbf{q}= [u_1,\ldots,u_d,p,T,(k,\omega),(\text{particle variables}),Y_1,\ldots,Y_{N-1}]^T.\] The physical carrier density \(\rho\) is constant. The stored conservative-like carrier components are \[\begin{aligned} W_{u_i} &= \rho u_i,\\ W_p &= p,\\ W_T &= \rho H-p =\rho\left(h(T,\mathbf{Y})+\frac12\lVert\mathbf{u}\rVert^2+k\right),\\ W_k &= \rho k,\qquad W_\omega=\rho\omega,\\ W_{Y_s}&=\rho Y_s. \end{aligned}\](4.1) The \(k\) contribution in \(W_T\) is present when GE-RANS is active.

Equation (4.1) is an important implementation detail: the pressure slot stores \(p\). It is not \(p/c_{ac}^{2}\). The artificial-compressibility scale enters the pseudo-time pressure equation described next.

4.2 Artificial pressure equation

The steady incompressible constraint is \[\nabla\cdot(\rho\mathbf{u})=0.\] CMPS creates a pseudo-time pressure-velocity coupling through \[\frac{1}{c_{ac}^{2}}\frac{\partial p}{\partial\tau} +\nabla\cdot(\rho\mathbf{u})=0, \](4.2) where \(\tau\) is pseudo time. Consequently the pressure-row pseudo-time Jacobian contains \(1/c_{ac}^{2}\).

The local artificial speed is bounded away from zero, \[c_{ac}=\max\left(U_{\mathrm{ref,min}},\;r_{ac}\lVert\mathbf{u}\rVert\right), \](4.3) with a positive reference-speed floor. The implementation also incorporates boundary-pressure and physical dual-time information when preparing local preconditioning/time-scale data. The floor in Eq. (4.3) prevents the pseudo-acoustic scale from vanishing at a stagnation point.

4.3 Interior inviscid flux

Let \(m_f\) be the numerical carrier mass flux per unit area and \(p_f\) the numerical mechanical pressure. For an interior face with outward unit normal \(\mathbf{n}\), the current incompressible carrier flux has the structure \[\begin{aligned} F_{u_i} &= m_f u_{i,\mathrm{up}} + p_f n_i,\\ F_p &= m_f,\\ F_T &= m_f H_{\mathrm{up}},\\ F_{Y_s} &= m_fY_{s,\mathrm{up}}. \end{aligned}\](4.4) GE-RANS uses the same numerical mass flux for \(k\) and \(\omega\) advection. Its mechanical pressure includes the modeled \(2\rho k/3\) contribution and the carrier enthalpy contains \(k\).

Using one numerical mass flux in continuity, momentum advection, total enthalpy, turbulence, and species transport is important for discrete consistency. The upwind side for transported quantities is selected from the sign of \(m_f\).

4.4 Artificial-compressibility AUSM⁺-up splitting

The active incompressible carrier scheme follows the AUSM⁺-up construction (Liou 2006), but its characteristic speed is the artificial speed rather than the thermodynamic sound speed. A shared face value \[c_{ac,f}=\max(c_{ac,L},c_{ac,R})\] defines the normal pseudo-Mach numbers \[M_L=\frac{u_{n,L}}{c_{ac,f}}, \qquad M_R=\frac{u_{n,R}}{c_{ac,f}}.\] The face Mach number is assembled from fourth-order split functions plus the AUSM⁺-up pressure-difference correction, \[M_f=\mathcal{M}^{+}_{4}(M_L)+\mathcal{M}^{-}_{4}(M_R)+M_p,\] and the face pressure is assembled from fifth-order pressure splits plus the velocity-difference correction, \[p_f=\mathcal{P}^{+}_{5}(M_L)p_L +\mathcal{P}^{-}_{5}(M_R)p_R+p_u.\] The implementation uses the standard AUSM⁺-up constants associated with its active path (\(\beta=1/8\), \(K_p=0.25\), \(K_u=0.75\), and \(\sigma=1\)). The pressure splitting is gauge-centred: adding the same constant \(C\) to both pressure states leaves the numerical mass flux unchanged and produces \(p_f\mapsto p_f+C\).

For constant density, \[m_f=\rho\,c_{ac,f}\,M_f,\] with the sign of \(m_f\) controlling the upwind state in Eq. (4.4).

4.5 Gauge pressure

The primitive pressure is not clipped to a positive absolute-pressure floor in the incompressible update. Negative gauge pressure is therefore representable. For a uniform gauge shift \[p\rightarrow p+C,\] the face mass flux remains invariant and the mechanical face pressure shifts by the same constant. This gives the expected gauge covariance in the momentum/continuity subsystem.

The energy variable deserves a qualification. Because \(H\) contains \(p/\rho\), an individual face enthalpy flux changes under a gauge shift. At a converged incompressible solution that change cancels with discrete continuity when the complete residual is considered. During finite pseudo-iterations, however, energy residual histories can depend on the selected pressure datum. Pressure reference and surface-force reference values should therefore be treated deliberately in gauge-pressure calculations.

4.6 Closed-domain pressure datum

A fully closed incompressible domain has the familiar arbitrary pressure constant. The current pseudo-pressure diagonal makes the linear update nonsingular during pseudo-time iteration and preserves the initialized gauge, but the implementation does not impose a separate pressure pin or a global mean-zero constraint. Consequently two independently initialized closed-domain runs need not converge to exactly the same additive pressure datum even when the velocity field is equivalent.

4.7 Physical-time treatment

For transient incompressible calculations, momentum, thermal, turbulence, particle, and species rows receive the selected BDF physical-time residual. The artificial pressure row does not receive a physical pressure-storage term. It remains the pseudo-time constraint used to converge each physical time step. The resulting algorithm is a dual-time method: physical accuracy is set by the BDF discretization, while artificial compressibility controls convergence of the inner pressure-velocity iterations.

4.8 Current implementation status

The present source path contains the dedicated incompressible selection, interior AUSM⁺-up flux, AD assembly, steady pseudo-time operation, dual-time terms, MPI worker execution, and agglomeration-multigrid dispatch. Formula-level tests of the audited flux reproduced uniform-flow consistency, face-orientation antisymmetry, constant-pressure shear, and gauge-shift invariance of mass flux to floating-point roundoff.

These observations are source/formula verification, not full runtime validation. The implementation review also identified open issues outside the core face formula, including validation bypass paths, unsupported TRANSMISSIVE assembly, restart history, moving-frame energy work, moving-wall work, particle drag-energy conservation, and non-adiabatic wall/CHT heat-balance problems. Chapter 20 lists these separately so that the mathematical formulation is not confused with release readiness.