20. Double-Flux Method for a Pressure-Based Segregated Solver#

This appendix describes the double-flux method used to reduce spurious pressure oscillations in multicomponent and real-fluid compressible flows. The main reference for the present formulation is the real-fluid double-flux treatment of Ma, Lv, and Ihme [MLI17], which extends the Abgrall-Karni and Billet-Abgrall double-flux idea to transcritical equations of state [AK01, BA03].

Ma et al. present the method for an explicit, density-based, conservative finite-volume algorithm. The solver described in this guide is pressure based and segregated: pressure is solved through a pressure-correction equation, velocity and fluxes are corrected, and total energy or total enthalpy is advanced in a separate scalar solve. The purpose of this appendix is to separate the thermodynamic idea from the density-based algorithm in which it is often introduced, then state the corresponding pressure-based equations and assumptions.

In one sentence:

During one nonlinear update, replace the nonlinear real-fluid energy-pressure relation by a cell-local frozen calorically-perfect-gas surrogate in the parts of the discretization that couple pressure, density, and energy-like transport; after that update, return to the true equation of state.

The word “surrogate” is essential. Double flux is not a new physical equation of state. It is a temporary, locally frozen thermodynamic representation used to keep an isobaric material interface from generating pressure through the discrete energy-pressure conversion.

20.1. Source of Pressure Oscillations#

For a general multicomponent real fluid, pressure is a nonlinear function of density, internal energy, and composition,

(20.1.1)#\[p = p(\rho,e,Y).\]

Across a material interface, pressure and velocity may be uniform while density, temperature, composition, sound speed, and internal energy vary sharply. A fully conservative density-based update transports mass, momentum, total energy, and species. Pressure is then recovered from the nonlinear EOS. If the local relation between energy and pressure changes abruptly across the interface, a conservative update that is reasonable for the conserved variables can still recover a nonuniform pressure.

Ma et al. analyze this mechanism by writing a general real-fluid state in the locally calorically-perfect-gas-like form

(20.1.2)#\[e = \frac{p v}{\gamma_{\mathrm{DF}}-1} + e_{0,\mathrm{DF}}, \qquad v = \frac{1}{\rho}.\]

The effective heat-capacity ratio is defined from the true thermodynamic sound speed,

(20.1.3)#\[\gamma_{\mathrm{DF}} = \frac{\rho a_s^2}{p},\]

and the reference internal energy is chosen so that the surrogate matches the current true-EOS state exactly,

(20.1.4)#\[e_{0,\mathrm{DF}} = e - \frac{p}{\rho(\gamma_{\mathrm{DF}}-1)}.\]

For a calorically perfect gas, these quantities reduce to constants. For a transcritical real fluid, they can vary strongly with pressure, density, temperature, and composition. It is often useful to introduce

(20.1.5)#\[\Gamma_{\mathrm{DF}} \equiv \frac{1}{\gamma_{\mathrm{DF}}-1},\]

so that Eq. (20.1.2) becomes

(20.1.6)#\[e = \Gamma_{\mathrm{DF}}\frac{p}{\rho}+e_{0,\mathrm{DF}}.\]

The pressure-equilibrium error is tied to jumps in the two state-dependent quantities

(20.1.7)#\[\Gamma_{\mathrm{DF}}, \qquad e_{0,\mathrm{DF}}.\]

Near pseudo-boiling or liquid-like to gas-like transitions, these jumps can be much larger than in ideal-gas mixing. This is why adding dissipation to a fully conservative method may not remove the pressure oscillation unless the interface is resolved on an impractically fine mesh.

20.2. Contact-Interface Balance#

The pressure-equilibrium argument is clearest for a one-dimensional contact-like state with uniform pressure and velocity. A first-order upwind update of continuity and momentum preserves uniform velocity. The question is whether the energy update preserves uniform pressure.

Using Eq. (20.1.6), the pressure change after one update can be written schematically as

(20.2.1)#\[\rho\,\delta e_{0,\mathrm{DF}} + p\,\delta \Gamma_{\mathrm{DF}} + \Gamma_{\mathrm{DF}}\,\delta p = \text{transport of spatial jumps in } e_{0,\mathrm{DF}} \text{ and } \Gamma_{\mathrm{DF}}.\]

Here \(\delta(\cdot)\) denotes the change over the update. The exact stencil indexing depends on the finite-volume discretization, but the interpretation is robust: pressure changes are driven by temporal and spatial variation of \(\Gamma_{\mathrm{DF}}\) and \(e_{0,\mathrm{DF}}\).

The double-flux idea is to hold \(\gamma_{\mathrm{DF}}\) and \(e_{0,\mathrm{DF}}\) fixed for each receiving cell over the update. With these quantities frozen, the cell no longer samples different nonlinear energy-pressure tangents during the same discrete update, and the contact update gives

(20.2.2)#\[\delta p = 0.\]

This does not make the global real-fluid EOS linear. It only makes the local energy-pressure conversion linear, and only for the duration of the frozen update.

20.3. Density-Based Double-Flux Form#

The conservative inviscid state used by Ma et al. is

(20.3.1)#\[\mathbf U = \left[ \rho,\; \rho \vec u,\; \rho E,\; \rho Y_1,\ldots,\rho Y_{N_s} \right]^T,\]

where

(20.3.2)#\[E = e + \frac{1}{2}\|\vec u\|^2.\]

For a cell whose frozen pair is \((\gamma_{\mathrm{DF}},e_{0,\mathrm{DF}})\), the double-flux total-energy relation is

(20.3.3)#\[\rho E_{\mathrm{DF}} = \frac{p}{\gamma_{\mathrm{DF}}-1} + \rho e_{0,\mathrm{DF}} + \frac{1}{2}\rho\|\vec u\|^2.\]

During a density-based Runge-Kutta stage, this relation is used when a total energy is needed for reconstruction or flux evaluation. After the conservative variables have been updated, pressure is recovered from the frozen relation,

(20.3.4)#\[p = (\gamma_{\mathrm{DF}}-1) \left[ \rho E - \rho e_{0,\mathrm{DF}} - \frac{1}{2}\rho\|\vec u\|^2 \right].\]

At the end of the time step, the thermodynamic state is synchronized again with the true EOS, and total energy is recomputed from the true thermodynamics. This synchronization is part of the method, not a cosmetic correction.

20.4. Cellwise Freezing and Quasi-Conservation#

The method is called “double flux” because an interior face generally contributes two different energy-like fluxes. If cells \(P\) and \(N\) share a face \(f\), the contribution to the row of \(P\) uses the frozen pair of \(P\), while the contribution to the row of \(N\) uses the frozen pair of \(N\).

For total energy, the same face state is viewed through two cell-local closures:

(20.4.1)#\[(\rho E)_{f,\mathrm{DF}}^{P} = \frac{p_f}{\gamma_{\mathrm{DF},P}-1} + \rho_f e_{0,\mathrm{DF},P} + \frac{1}{2}\rho_f\|\vec u_f\|^2,\]

(20.4.2)#\[(\rho E)_{f,\mathrm{DF}}^{N} = \frac{p_f}{\gamma_{\mathrm{DF},N}-1} + \rho_f e_{0,\mathrm{DF},N} + \frac{1}{2}\rho_f\|\vec u_f\|^2.\]

In general,

(20.4.3)#\[(\rho E)_{f,\mathrm{DF}}^{P} \ne (\rho E)_{f,\mathrm{DF}}^{N}.\]

Mass, momentum, and species remain conservative in the finite-volume sense. Total energy is quasi-conservative because the energy-like flux is intentionally localized to the receiving cell’s frozen thermodynamic tangent. Ma et al. show that the resulting energy conservation error is tied to the same jumps in Eq. (20.1.7). For smooth solutions it decreases with mesh refinement. For shock-dominated problems, especially on coarse meshes, the error can affect shock speed and temperature. This is the main motivation for hybrid conservative / double-flux strategies [BJ21].

20.5. Translation to a Pressure-Based Solver#

In a density-based algorithm, pressure is recovered after a conservative update. In a pressure-based segregated algorithm, pressure is a primary unknown and density changes enter through the pressure-correction equation. Therefore, the translation is not to copy the explicit Runge-Kutta algorithm stage by stage. The translation is to preserve the same frozen thermodynamic relation inside the pressure-density-energy coupling of one nonlinear update.

Conceptually, a pressure-based double-flux update uses the following contract:

Start from a state that is consistent with the true EOS.
Evaluate the true sound speed and internal energy at that state.
Form one frozen pair \((\gamma_{\mathrm{DF}},e_{0,\mathrm{DF}})\) per cell using Eqs. (20.1.3) and (20.1.4).
Hold that pair fixed while assembling and solving the pressure-correction and energy-like scalar operators built from the same nonlinear reference state.
Use row-local double-flux energy-like face values for energy or enthalpy transport.
Use the frozen-surrogate pressure-density response in the pressure-correction terms that represent compressible density change.
After pressure, velocity, density, species, and energy-like variables have been updated, close the thermodynamic state again with the true EOS.

The frozen pair is cell local. It is not a face average, and it is not updated inside a linear solve. The true EOS remains the source of sound speed, internal energy, final thermodynamic properties, transport properties, and any chemistry or tabulated-thermodynamic closure.

20.6. Energy and Enthalpy Forms#

A pressure-based solver may advance either total energy or total enthalpy. Both are compatible with the same frozen surrogate. Define

(20.6.1)#\[\alpha_{\mathrm{DF}} \equiv \frac{\gamma_{\mathrm{DF}}}{\gamma_{\mathrm{DF}}-1}.\]

The frozen static enthalpy relation is

(20.6.2)#\[h_{\mathrm{DF}} = e_{0,\mathrm{DF}} + \alpha_{\mathrm{DF}}\frac{p}{\rho}.\]

The frozen total enthalpy relation is

(20.6.3)#\[H_{\mathrm{DF}} = e_{0,\mathrm{DF}} + \alpha_{\mathrm{DF}}\frac{p}{\rho} + \frac{1}{2}\|\vec u\|^2.\]

If density is needed from a frozen total-energy relation, the inversion is

(20.6.4)#\[\rho_{\mathrm{DF}} = \frac{p} {(\gamma_{\mathrm{DF}}-1) \left(E-\frac{1}{2}\|\vec u\|^2-e_{0,\mathrm{DF}}\right)}.\]

If density is needed from a frozen total-enthalpy relation, the inversion is

(20.6.5)#\[\rho_{\mathrm{DF}} = \alpha_{\mathrm{DF}} \frac{p} {H-\frac{1}{2}\|\vec u\|^2-e_{0,\mathrm{DF}}}.\]

These inversions are not new physical EOS calls. They are operator-local closures implied by the frozen double-flux relation. The final state is still closed by the true EOS.

20.7. Pressure-Correction Density Response#

Without double flux, the pressure-correction derivation uses the true acoustic density response

(20.7.1)#\[\rho' \approx \frac{p'}{a_s^2}.\]

This is the EOS-compatible response along an acoustic, entropy-fixed tangent of the true thermodynamic state.

With double flux active, the pressure-density relation used inside the frozen update is different. Holding the energy-like variable, velocity, \(\gamma_{\mathrm{DF}}\), and \(e_{0,\mathrm{DF}}\) fixed in Eq. (20.6.4) or Eq. (20.6.5) gives

(20.7.2)#\[\rho' \approx \frac{\rho}{p}\,p'.\]

Using Eq. (20.1.3), this can be written as

(20.7.3)#\[\frac{\rho}{p} = \frac{\gamma_{\mathrm{DF}}}{a_s^2}.\]

Thus the double-flux pressure-correction density response is the frozen-surrogate slope \(\rho/p\), equivalently \(\gamma_{\mathrm{DF}}/a_s^2\), not the pure acoustic slope \(1/a_s^2\). This is consistent with the pressure recovery relation in the density-based double-flux method: during the frozen update, pressure and density are related through the local calorically-perfect-gas surrogate.

This distinction does not abandon the true EOS. The sound speed used to define \(\gamma_{\mathrm{DF}}\) comes from the true EOS at the reference state, and the final thermodynamic state is again closed with the true EOS.

20.8. Face Fluxes in the Pressure-Based Form#

The pressure-based solver keeps a single conservative mass flux at a face. The double-flux modification is applied to the energy-like quantity carried by that mass flux.

For a total-enthalpy form, the row-local face value for cell \(P\) is

(20.8.1)#\[H_{f,\mathrm{DF}}^{P} = e_{0,\mathrm{DF},P} + \alpha_{\mathrm{DF},P} \frac{p_f}{\rho_f} + \frac{1}{2}\|\vec u_f\|^2.\]

The row-local energy-like convective contribution is then

(20.8.2)#\[\dot m_f\,H_{f,\mathrm{DF}}^{P}.\]

The neighboring row uses the neighboring frozen pair,

(20.8.3)#\[H_{f,\mathrm{DF}}^{N} = e_{0,\mathrm{DF},N} + \alpha_{\mathrm{DF},N} \frac{p_f}{\rho_f} + \frac{1}{2}\|\vec u_f\|^2.\]

The pressure, density, velocity, and composition at the face should be reconstructed consistently with the inviscid flux used by the solver. The double-flux change is the thermodynamic conversion from that face state to the transported energy-like quantity. It is not a separate momentum flux, not a separate species flux, and not a replacement for the true EOS used to define the physical state before and after the update.

20.9. Assumptions and Limits#

The contact-pressure argument assumes that the problematic part of the update is the nonlinear conversion among pressure, density, and energy at a material interface. Composition may change by species transport and reactions, but the local thermodynamic conversion is frozen over the update being assembled.

The frozen surrogate requires nondegenerate states:

(20.9.1)#\[p > 0, \qquad \rho > 0, \qquad \gamma_{\mathrm{DF}} \ne 1,\]

and, when density is recovered from energy-like variables,

(20.9.2)#\[E-\frac{1}{2}\|\vec u\|^2-e_{0,\mathrm{DF}} > 0, \qquad H-\frac{1}{2}\|\vec u\|^2-e_{0,\mathrm{DF}} > 0.\]

If these conditions fail, the state is outside the useful range of the local frozen surrogate. A robust method should limit, reject, or locally fall back to a conservative / true-EOS treatment rather than continuing with an invalid density.

Double flux also does not provide positivity preservation, bounded species, entropy stability, or shock robustness by itself. Ma et al. combine the method with an entropy-stable hybrid correction for realizability and minimum-entropy control. In that setting, the entropy measure is built from the effective heat-capacity ratio because each cell is temporarily represented by a locally frozen calorically-perfect-gas surrogate.

20.10. Hybrid Conservative / Double-Flux Strategies#

The double-flux method is most useful where pressure equilibrium at a material or transcritical interface is more important than strict total-energy conservation over one under-resolved update. It is less attractive in shock-dominated regions where exact energy conservation is important for shock speed and post-shock temperature.

A common hybrid idea is to blend or switch between a conservative energy flux and a double-flux energy-like flux:

(20.10.1)#\[F_E = \theta_f F_E^{\mathrm{conservative}} + (1-\theta_f)F_E^{\mathrm{DF}}, \qquad 0 \le \theta_f \le 1.\]

Here \(\theta_f\) should be close to zero where thermodynamic-interface pressure oscillations dominate, and close to one where strict conservative shock behavior is more important. Boyd and Jarrahbashi [BJ21] discuss hybrid strategies based on jumps in \(\Gamma_{\mathrm{DF}}\) and \(e_{0,\mathrm{DF}}\), the same quantities identified in Eq. (20.1.7). A thermodynamic jump sensor identifies where double flux may be needed, but a shock-aware sensor or blending rule is usually needed as well.

20.11. Summary#

The double-flux method freezes a local calorically-perfect-gas surrogate,

(20.11.1)#\[p=(\gamma_{\mathrm{DF}}-1)\rho(e-e_{0,\mathrm{DF}}),\]

with

(20.11.2)#\[\gamma_{\mathrm{DF}}=\frac{\rho a_s^2}{p}, \qquad e_{0,\mathrm{DF}}= e-\frac{p}{\rho(\gamma_{\mathrm{DF}}-1)}.\]

The frozen pair removes the local temporal and spatial variation in the energy-pressure tangent that generates pressure at an isobaric material interface. The method is quasi-conservative in total energy because neighboring cells can use different frozen energy-like fluxes at the same face.

In a pressure-based segregated solver, the same thermodynamic idea appears as a frozen pressure-density-energy closure over one nonlinear update. The pressure-correction density response uses the frozen-surrogate slope \(\rho/p=\gamma_{\mathrm{DF}}/a_s^2\), energy or enthalpy transport uses row-local double-flux face values, and the updated thermodynamic state is then closed again with the true EOS.

20.12. References#

[AK01]

R. Abgrall and S. Karni. Computations of compressible multifluids. Journal of Computational Physics, 169(2):594–623, 2001. doi:10.1006/jcph.2000.6685.

[BA03]

G. Billet and R. Abgrall. An adaptive shock-capturing algorithm for solving unsteady reactive flows. Computers and Fluids, 32(10):1473–1495, 2003.

[BJ21] (1,2)

B. Boyd and D. Jarrahbashi. A diffuse-interface method for reducing spurious pressure oscillations in multicomponent transcritical flow simulations. Computers and Fluids, 222:104924, 2021. doi:10.1016/j.compfluid.2021.104924.

[MLI17]

P.C. Ma, Y. Lv, and M. Ihme. An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows. Journal of Computational Physics, 340:330–357, 2017. doi:10.1016/j.jcp.2017.03.022.