24. Pressure-Correction Equation Forms

Chapter 6 derives the pressure-correction row from the continuity defect and the face-mass-flux correction. The present appendix does not repeat that derivation. Instead, it covers the different ways that the pressure-correction equation can be expressed and how those changes manifest in the resulting algebraic equation.

This appendix compares alternative forms of the same pressure-correction row and makes the coefficient bookkeeping explicit. For the original derived form of the pressure-correction equation, see Pressure-Correction Equation on Unstructured Grids. The discussion below is intentionally written in generic notation first.

Four practical points are worth keeping in view from the outset:

1. all families still target one cellwise algebraic equation for \(p^{'}\)

2. the family choice mainly changes things such as the face coefficient multiplying \(p_N^{'}-p_P^{'}\) or the row normalization

3. compressible families may add an upwind acoustic density-correction term

4. the species-summed volumetric form is special because its normalization enters through species-weighted sums rather than by dividing the entire mixture row by \(\rho\)

The appendix follows the same reading pattern used in Pressure-Correction Equation on Unstructured Grids:

1. restate the common pressure-correction row,

2. isolate the face coefficient that comes from the pressure-difference flux correction,

3. show how the face-density normalization, the rowwise density normalization, and SIMPLE/SIMPLEC modify that coefficient,

4. collect the final diagonal, neighbor, and right-hand-side coefficients for each implemented family.

24.1. Reference Equation Form

The common assembled pressure-correction row is

(24.1.1)\[a_P p_P^{'} = \sum_{nbr} a_{nbr} p_nbr^{'} + b_P\]

This is the same positive-neighbor form used in (6.5.6). The form used when constructing the linear solve system is the sign-absorbed version of the same equation:

(24.1.2)\[{}^{(p')}\!D_P p_P^{'} + \sum_{nbr} {}^{(p')}\!A_{nbr} p_nbr^{'} = {}^{(p')}\!B_P, \qquad {}^{(p')}\!A_{nbr} = -a_{nbr}.\]

The discussion below stays primarily in the assembled form of Eq. (24.1.1), because that is the clearest place to compare the alternative pressure-correction forms. The equivalent sign-absorbed storage convention is the same one summarized in Pressure-Correction Matrix Storage.

The Main Question This Appendix Answers

For one interior face between cells \(P\) and \(N\), the practical question is: what coefficient finally multiplies the pressure-difference term in the \(P\) row, and how does that coefficient change when the formulation switches forms?

Everything that follows is answering that one question in the same style used in Pressure-Correction Equation on Unstructured Grids.

24.2. Common Algebraic Backbone

Begin from the same cellwise continuity balance used in (6.1.1):

(24.2.1)\[\left[ \frac{\rho_P^{k+1}-\rho_P^{n}}{\Delta t} \right]\Omega_P + \sum_{f\in\partial P} \left( \hat{m}_{f}^{*} + \dot{m}_{f}^{'} \right) =0.\]

Here \(\hat{m}_{f}^{*}\) is the known provisional face flux written in the orientation of the row for cell \(P\), and \(\dot{m}_{f}^{'}\) is the unknown pressure-induced correction.

As in Pressure-Correction Equation on Unstructured Grids, the compressible families use the acoustic EOS linearization

(24.2.2)\[\rho_P^{'} \approx \frac{p_P^{'}}{a_P^2},\]

Writing

(24.2.3)\[\rho_P^{k+1} = \rho_P^k + \rho_P^{'},\]

and substituting Eq. (24.2.2) into the unsteady term gives

(24.2.4)\[\left[ \frac{\rho_P^{k+1}-\rho_P^{n}}{\Delta t} \right]\Omega_P \approx \frac{\Omega_P}{\Delta t\,a_P^2}p_P^{'} + \frac{\rho_P^k-\rho_P^n}{\Delta t}\Omega_P.\]

The first term is proportional to the unknown pressure correction, so it becomes the implicit unsteady coefficient in the diagonal. The second term depends only on already known iteration states, so it is treated as a lagged source contribution on the right-hand side.

At the face level, all of the standard mixture-based forms can be viewed as using the same base pressure-difference correction pattern,

(24.2.5)\[\dot{m}_{f}^{'} = - d_f\left(p_N^{'}-p_P^{'}\right) + \text{optional compressible additions}\]

where \(d_f\) denotes the abstract base face coefficient introduced in Eq. (6.4.2).

To isolate the part common to all forms, first ignore the extra compressible density-correction terms that are added later for the compressible families. With only that base pressure-difference term retained, one interior face contributes

(24.2.6)\[\dot{m}_{f}^{'} = d_f p_P^{'} - d_f p_N^{'}\]

to the \(P\) row. This is exactly the same diagonal/neighbor bookkeeping step emphasized in Pressure-Correction Equation on Unstructured Grids.

Collecting those face contributions for all the faces of cell \(P\) gives the same base algebraic pattern as (6.5.6):

(24.2.7)\[a_P^{(0)} p_P^{'} = \sum_{nbr} a_{nbr}^{(0)} p_{nbr}^{'} + b_P^{(0)}\]

with

(24.2.8)\[a_{nbr}^{(0)} = d_f\]

(24.2.9)\[a_P^{(0)} = \sum_{f} d_f + \sum_{bnd} d_b + {}^{(p')}\!a_{P,\mathrm{time}},\]

and

(24.2.10)\[b_P^{(0)} = -\sum_{f} \hat{m}_{f}^{*} + b_{P,\mathrm{time}} + b_{P,\mathrm{bnd}}.\]

So the algebraic spine is the same as in Pressure-Correction Equation on Unstructured Grids. The remaining question is how the coefficient \(d_f\) from that derivation is rescaled into the family-dependent coefficient actually used by each pressure-correction form.

24.3. The Face-Coefficient Pipeline

Up to this point, \(d_f\) has been the abstract base coefficient multiplying the pressure-difference term in the face-flux correction. The purpose of the present section is to unpack the successive algebraic stages that turn that base coefficient into the form-specific coefficient used in the assembled row.

To keep the notation systematic, this appendix uses the following convention:

• \(d_f\) denotes the base pressure-difference face coefficient from Eq. (6.4.2),

• \(C\) denotes the same coefficient after face-density normalization but before SIMPLE/SIMPLEC scaling,

• \(D\) denotes the same coefficient after SIMPLE/SIMPLEC scaling,

• the superscript \(\mathrm{row}\) means “the coefficient in the normalized row of cell \(P\) after rowwise normalization.”

With that convention, one interior-face coefficient is related by the chain

(24.3.1)\[d_f \;\longrightarrow\; C_f = \frac{d_f}{\hat{\rho}^{(p)}_f} \;\longrightarrow\; D_f = \psi_f C_f \;\longrightarrow\; D_{P,f}^{(\mathrm{row})} = \frac{D_f}{\hat{\rho}^{(n)}_P} = \frac{\psi_f d_f}{\hat{\rho}^{(p)}_f\,\hat{\rho}^{(n)}_P}.\]

So \(d_f\), \(C_f\), \(D_f\), and \(D_{P,f}^{(\mathrm{row})}\) are not different physical effects. They are four algebraic versions of the same face coefficient viewed before normalization, after face normalization, after SIMPLE/SIMPLEC scaling, and after rowwise normalization in the row of cell \(P\).

24.3.1. Form-Control Switches

To write several mixture-based pressure-correction forms with one compact algebraic template, introduce two 0/1 form-control switches:

(24.3.1.1)\[\hat{\rho}^{(p)}_f = (1-\sigma_p)\rho_f + \sigma_p, \qquad \hat{\rho}^{(n)}_P = (1-\sigma_\rho)\rho_P + \sigma_\rho,\]

with

(24.3.1.2)\[\sigma_p \in\{0,1\}, \qquad \sigma_\rho \in\{0,1\}.\]

These switches are bookkeeping devices, not new physical variables. The switch \(\sigma_p\) controls whether the face density remains inside the pressure-difference coefficient, and \(\sigma_\rho\) controls whether the assembled row is kept in mixture-mass form or divided by the cell density. The concrete option names that realize these choices are summarized later in the code-mapping section.

Because each switch takes only the values 0 or 1, the pair \((\sigma_p,\sigma_\rho)\) produces four algebraic combinations. To keep the names systematic, this appendix uses a two-part label: the first part describes what happens to the face coefficient, and the second part describes what happens to the assembled row.

\(\sigma_p\)	\(\sigma_\rho\)	Systematic name	Meaning
1	1	retained-face mass-row form	keep the face density in the pressure-difference coefficient, and keep the assembled row in mass form
0	1	divided-face mass-row form	divide the face density out of the pressure-difference coefficient, and keep the assembled row in mass form
1	0	retained-face volumetric-row form	keep the face density in the pressure-difference coefficient, and divide the assembled row by the cell density
0	0	divided-face volumetric-row form	divide the face density out of the pressure-difference coefficient, and divide the assembled row by the cell density

The first three of these combinations appear as named families in the present solver. The fourth is algebraically well defined, but it is not one of the named implementation paths documented in this appendix.

24.3.2. Face-Density Normalization

The face-flux derivation in Pressure-Correction Equation on Unstructured Grids introduced the abstract coefficient \(d_f\) through the prototype face-flux correction

(24.3.2.1)\[\dot{m}_f^{'} = -d_f\left(p_N^{'}-p_P^{'}\right),\]

in Eq. (6.4.2). That same \(d_f\) is the starting coefficient here; the present appendix simply tracks how the selected form rescales it.

With that notation, the face-density-normalized coefficient is

(24.3.2.2)\[C_f = \frac{d_f}{\hat{\rho}^{(p)}_f}.\]

For the standard mixture-based forms, that base coefficient factors as

(24.3.2.3)\[d_f = \rho_f\,\widetilde{C}_f,\]

where \(\widetilde{C}_f\) means “everything in the base momentum-interpolation coefficient except the explicit face-density factor.” In other words, \(\widetilde{C}_f\) collects the remaining geometric, momentum-response, and variant-specific consistency factors once \(\rho_f\) has been factored out.

So Eq. (24.3.2.2) becomes

(24.3.2.4)\[C_f = \frac{\rho_f\,\widetilde{C}_f}{(1-\sigma_p)\rho_f+\sigma_p}.\]

Equivalently, the same result can be written in case form as

(24.3.2.5)\[\begin{split}C_f = \begin{cases} \rho_f\widetilde{C}_f, & \sigma_p=1, \\ \widetilde{C}_f, & \sigma_p=0 . \end{cases}\end{split}\]

• \(\sigma_p=1\): this is the retained-face branch, so the pressure-correction coefficient keeps the factor \(\rho_f\).

• \(\sigma_p=0\): this is the divided-face branch, so the pressure-correction coefficient uses only \(\widetilde{C}_f\).

So \(\sigma_p\) changes only one thing: whether the pressure-difference term is weighted by the face density before any rowwise normalization is applied.

24.3.3. SIMPLE and SIMPLEC Scaling

After the face-density normalization choice has been made, the momentum-correction part is scaled by

(24.3.3.1)\[\psi_f = 1+\kappa_f, \qquad \kappa_f = \frac{\alpha_u}{1-\alpha_u+\delta_f},\]

while SIMPLE uses \(\psi_f=1\).

So the pressure-difference coefficient before any rowwise density normalization is

(24.3.3.2)\[D_f = \psi_f C_f.\]

24.3.4. What Rowwise Normalization Means

Before specializing to density-based normalization, it helps to state the general algebraic idea. If one assembled row is written as

(24.3.4.1)\[a_P x_P = \sum_{N\in\mathcal{N}(P)} a_{PN} x_N + b_P,\]

then rowwise normalization by any nonzero cellwise factor \(s_P\) means dividing the entire row by that same factor:

(24.3.4.2)\[\frac{a_P}{s_P} x_P = \sum_{N\in\mathcal{N}(P)} \frac{a_{PN}}{s_P} x_N + \frac{b_P}{s_P}.\]

So the diagonal, every neighbor coefficient, and the right-hand side are all scaled together. The equation represented by the row is unchanged; what changes is the normalization and therefore the interpretation and conditioning of the assembled coefficients.

24.3.5. Rowwise Density Normalization

The switch \(\sigma_\rho\) acts at the row level rather than at the base face-coefficient level. In other words, it divides the entire assembled row, not just the pressure-difference term. If \(a_{P\star}^{(\text{mass})}\) and \(b_P^{(\text{mass})}\) denote the coefficients before this row normalization, then

(24.3.5.1)\[a_{P\star}^{(\mathrm{row})}=\frac{a_{P\star}^{(\text{mass})}}{\hat{\rho}^{(n)}_P}, \qquad b_P^{(\mathrm{row})}=\frac{b_P^{(\text{mass})}}{\hat{\rho}^{(n)}_P}.\]

Therefore:

• \(\sigma_\rho=1\): \(\hat{\rho}^{(n)}_P=1\), so the mixture mass-form row is retained.

• \(\sigma_\rho=0\): the entire row is divided by \(\rho_P\), giving the normalized mixture volumetric form.

This is the continuous distinction between

(24.3.5.2)\[\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0\]

and

(24.3.5.3)\[\frac{1}{\rho}\left[\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})\right]=0.\]

The key structural point is that \(\sigma_p\) and \(\sigma_\rho\) do not do the same job. The switch \(\sigma_p\) changes only the pressure-difference face coefficient. The switch \(\sigma_\rho\) scales the entire row, including the right-hand side and any unsteady term.

Final Face Coefficient Used in Row \(P\)

Putting the four steps together, the coefficient that actually multiplies the pressure-difference term in the row of cell \(P\) is

(24.3.5.4)\[D_{P,f}^{(\mathrm{row})} = \frac{\psi_f d_f}{\hat{\rho}^{(p)}_f\,\hat{\rho}^{(n)}_P}.\]

For the standard base form \(d_f=\rho_f\widetilde{C}_f\), this becomes

(24.3.5.5)\[D_{P,f}^{(\mathrm{row})} = \psi_f \frac{\rho_f\widetilde{C}_f} {\left[(1-\sigma_p)\rho_f+\sigma_p\right]\left[(1-\sigma_\rho)\rho_P+\sigma_\rho\right]}.\]

This is the compact algebraic answer to the question “what final face coefficient does this pressure-correction form use in the \(P\) row?”

24.3.6. Four Switch Combinations

When the family uses the standard base face coefficient \(\rho_f\widetilde{C}_f\), Eq. (24.3.5.5) immediately gives the four algebraic switch combinations:

(24.3.6.1)\[\begin{split}D_{P,f}^{(\mathrm{row})} = \begin{cases} \psi_f\,\rho_f\widetilde{C}_f, & \text{retained-face mass-row form}, \\ \psi_f\,\widetilde{C}_f, & \text{divided-face mass-row form}, \\ \psi_f\,\dfrac{\rho_f\widetilde{C}_f}{\rho_P}, & \text{retained-face volumetric-row form}, \\ \psi_f\,\dfrac{\widetilde{C}_f}{\rho_P}, & \text{divided-face volumetric-row form}. \end{cases}\end{split}\]

The present solver gives distinct family names only to the first three combinations. Depending on flow regime, the retained-face mass-row form appears as the constant-density mass form or the mixture mass form. The divided-face mass-row form is the face-density-divided form, and the retained-face volumetric-row form is the mixture volumetric form. The divided-face volumetric-row form is included here to complete the four-way switch map, but it is not one of the named implementation paths discussed later in this appendix.

The first two lines differ only by the removal of the face-density factor. So, in the constant-density incompressible limit, the divided-face mass-row form is exactly the usual divergence-form version of the same pressure-correction equation. The third line keeps the face-density factor in the pressure-difference term but divides the entire row by the cell density, which is why it represents a normalized volumetric mixture-continuity form.

24.4. Form Families

At the theory level, the supported named pressure-correction families can be summarized as follows:

Flow regime	Named family	\(\sigma_p\)	\(\sigma_\rho\)	Continuity / row form represented
incompressible	constant-density mass form	1	1	constant-density mass-form row
incompressible	face-density-divided form	0	1	same row written without the face-density factor; in the constant-density limit this is the usual divergence form, \(\nabla\cdot \mathbf{u}=0\)
incompressible	compressible bookkeeping path	1	1	mass-form bookkeeping only; finite-sound-speed additions remain inactive because the flow mode is incompressible
compressible	mixture mass form	1	1	mixture mass continuity, \(\partial_t\rho+\nabla\cdot(\rho\mathbf{u})=0\), with acoustic additions
compressible	mixture volumetric form	1	0	mixture row divided by \(\rho_P\), i.e. the normalized volumetric form of Eq. (24.3.5.3)
compressible	species-summed volumetric form	1	1	normalized sum of the species equations; the normalization is carried by the species-weighted sums rather than by dividing the whole mixture row by \(\rho\)

The first five rows fit directly into the four-switch map introduced above. The species-summed volumetric form is the special case discussed later, where the volumetric normalization is carried by species-weighted factors rather than by the switch \(\sigma_\rho\). The closing code-mapping section lists the concrete option names used by the implementation.

24.5. Additional Compressible Density-Correction Term at Faces

The base pressure-difference term described above is present in both incompressible and compressible families. Compressible families add one more ingredient: the density correction carried by the provisional relative face flux.

A compact magnitude estimate for that extra contribution is

(24.5.1)\[K_f = \chi_c\, \frac{\left|\dot{V}_{f,\mathrm{rel}}^{*}\right|}{a_{\text{up}}^2},\]

where \(a_{\text{up}}\) is the upwind thermodynamic sound speed returned by the EOS and \(\chi_c\) is the coefficient multiplying the implicit compressible density correction in the selected form.

The underlying density linearization is the same acoustic Taylor expansion used in Pressure-Correction Equation on Unstructured Grids. If

(24.5.2)\[\widetilde{\rho} \equiv \rho(p^*+p',s^*,Y^*), \qquad \rho^* \equiv \rho(p^*,s^*,Y^*),\]

then

(24.5.3)\[\widetilde{\rho} = \rho^* + \left.\left( \frac{\partial \rho}{\partial p} \right)_{s,Y}\right|_{(p^*,\,s^*,\,Y^*)} p' + O\!\left((p')^2\right),\]

with the acoustic identity

(24.5.4)\[\left.\left( \frac{\partial \rho}{\partial p} \right)_{s,Y}\right|_{(p^*,\,s^*,\,Y^*)} = \frac{1}{a^2}.\]

Therefore the density correction is

(24.5.5)\[\rho' \equiv \widetilde{\rho}-\rho^* \approx \left.\left( \frac{\partial \rho}{\partial p} \right)_{s,Y}\right|_{(p^*,\,s^*,\,Y^*)} p' = \frac{p'}{a^2}.\]

The appearance of \(s\) here does not mean that the pressure-correction equation is being written in an entropy formulation. It only identifies the acoustic tangent of the EOS used in the linearization. This is an EOS-compatible pressure-correction linearization, not a hardwired ideal-gas \((\partial\rho/\partial p)_{T,Y}\) assumption.

Assumptions Behind The Compressible Additions

Three assumptions are built into the compressible pressure-correction term:

1. only the pressure-sensitive part of the density change is treated implicitly in the \(p^{'}\) solve;

2. the relevant thermodynamic slope is the acoustic derivative \(1/a^2\);

3. the extra density-correction term is upwinded with the provisional relative face flux.

So this is a pressure-correction linearization step, not a claim that all density variation in the full nonlinear problem is exhausted by \(p^{'}\).

For a face oriented from cell \(P\) to cell \(N\), define the sign-split upwind coefficients

(24.5.6)\[K_{P,f}= \chi_c \frac{\max\left(\dot{V}^{*}_{f,\mathrm{rel}},0\right)}{a_P^2}, \qquad K_{N,f}= \chi_c \frac{\max\left(-\dot{V}^{*}_{f,\mathrm{rel}},0\right)}{a_N^2},\]

where \(\dot{V}^{*}_{f,\mathrm{rel}}\) is the provisional relative volumetric flux. The face mass-flux correction is then written as

(24.5.7)\[\dot{m}_f^{'} = -D_f\left(p_N^{'}-p_P^{'}\right) + K_{P,f}p_P^{'} - K_{N,f}p_N^{'}.\]

Expanding this in the row of cell \(P\) gives the final facewise coefficients

(24.5.8)\[a_{P,f}^{(\mathrm{comp})} = \frac{D_f+K_{P,f}}{\hat{\rho}^{(n)}_P}, \qquad a_{PN,f}^{(\mathrm{comp})} = \frac{D_f+K_{N,f}}{\hat{\rho}^{(n)}_P}.\]

So the compressible additions do not create a different matrix structure. They simply augment the diagonal of the upwind row and the neighbor coefficient of the downwind row, exactly as described in (6.6.5).

The unsteady part is normalized in the same rowwise way. For example, the explicit BDF-type defect term has the form

(24.5.9)\[b_{P,\text{unsteady}} \propto \frac{\rho_P^n\Omega_P^n-\rho_P^k\Omega_P^k}{\Delta t\,\hat{\rho}^{(n)}_P}.\]

For ideal gases, the acoustic derivative reduces to

(24.5.10)\[\frac{1}{a^2}=\frac{1}{\gamma R T}=\frac{\rho}{\gamma p}, \qquad \frac{1}{\rho a^2}=\frac{1}{\gamma p},\]

so the older ideal-gas formulas are recovered as special cases of the EOS-generic acoustic form above.

24.6. Species-Summed Volumetric Form

The third family is a multicomponent volumetric form derived from a normalized sum of the species equations rather than from mixture continuity. It does not obtain its normalization by dividing the entire row by the mixture density. Instead, the normalization is carried by species-weighted cell and face sums that appear naturally in the derivation.

The same pattern used in Pressure-Correction Equation on Unstructured Grids is used again: start from the governing balance, split the starred and corrected terms, define the cell and face factors, and then collect the final row.

24.6.1. Start from the Normalized Species Balance

For each species \(i\), write the transport equation as

(24.6.1.1)\[\frac{\partial(\rho Y_i)}{\partial t} + \nabla\cdot(\rho \mathbf{u} Y_i) = \mathcal{D}_i + \mathcal{S}_i,\]

where \(\mathcal{D}_i\) denotes the diffusive part of the species equation and \(\mathcal{S}_i\) denotes chemistry or other production/source terms.

Dividing by the local pure-species density \(\rho_i=1/v_i\) and summing over all species gives

(24.6.1.2)\[\sum_{i=1}^{N_s} \frac{1}{\rho_i} \left[ \frac{\partial(\rho Y_i)}{\partial t} + \nabla\cdot(\rho \mathbf{u} Y_i) \right] = \sum_{i=1}^{N_s}\frac{\mathcal{D}_i+\mathcal{S}_i}{\rho_i}.\]

Integrating over a control volume \(P\) and writing the convective term with the face mass flux \(\dot{m}_{f}\), understood in the orientation of the row for cell \(P\), gives

(24.6.1.3)\[\frac{\Omega_P}{\Delta t} \sum_{i=1}^{N_s} \frac{(\rho Y_i)_P-(\rho Y_i)_P^n}{\rho_{i,P}} + \sum_{f\in\partial P} \dot{m}_{f} \sum_{i=1}^{N_s}\frac{Y_{i,f}^{\mathrm{up}}}{\rho_{i,P}} = b_{P,\mathrm{sp}}^{\mathrm{lag}},\]

with

(24.6.1.4)\[b_{P,\mathrm{sp}}^{\mathrm{lag}} \equiv \Omega_P \sum_{i=1}^{N_s} \frac{\mathcal{D}_{i,P}^{\mathrm{lag}}+\mathcal{S}_{i,P}^{\mathrm{lag}}}{\rho_{i,P}}.\]

Here \(Y_{i,f}^{\mathrm{up}}\) denotes the face composition used by the active species-flux reconstruction. The rowwise detail to keep in mind is that the \(P\) row uses the pure-species densities from cell \(P\), while the neighboring row uses the same face composition with the pure-species densities from cell \(N\).

24.6.2. Introduce the Pressure-Correction Split

As in Pressure-Correction Equation on Unstructured Grids, write

(24.6.2.1)\[\mathbf{u}=\mathbf{u}^{*}+\mathbf{u}^{'}, \qquad p=p^{*}+p^{'}, \qquad \rho=\rho^{*}+\rho^{'}, \qquad Y_i^{'}\approx 0.\]

Substituting Eq. (24.6.2.1) into Eq. (24.6.1.3) gives

(24.6.2.2)\[\frac{\Omega_P}{\Delta t} \sum_{i=1}^{N_s} \frac{(\rho^{*}Y_i)_P-(\rho^n Y_i^n)_P}{\rho_{i,P}} + \frac{\Omega_P}{\Delta t} \sum_{i=1}^{N_s} \frac{\rho_P^{'}Y_{i,P}}{\rho_{i,P}} + \sum_{f\in\partial P}\dot{m}_{f}^{*} \sum_{i=1}^{N_s}\frac{Y_{i,f}^{\mathrm{up}}}{\rho_{i,P}} + \sum_{f\in\partial P}\dot{m}_{f}^{'} \sum_{i=1}^{N_s}\frac{Y_{i,f}^{\mathrm{up}}}{\rho_{i,P}} = b_{P,\mathrm{sp}}^{\mathrm{lag}}.\]

The structure is now exactly parallel to Pressure-Correction Equation on Unstructured Grids: an explicit lagged unsteady term, an implicit density-correction term, the known starred face-flux defect, and the corrected face-flux contribution.

24.6.3. Species-Weighted Normalization Factors

At the cell center,

(24.6.3.1)\[F_P=\sum_{i=1}^{N_s} Y_{i,P} v_{i,P}.\]

At an interior face between cells \(P\) and \(N\),

(24.6.3.2)\[F_{f,P}=\sum_{i=1}^{N_s} Y_{i,f}^{\mathrm{up}} v_{i,P}, \qquad F_{f,N}=\sum_{i=1}^{N_s} Y_{i,f}^{\mathrm{up}} v_{i,N}.\]

At a boundary face, this reduces to the one-sided factor \(F_{b,P}\).

With these definitions, Eq. (24.6.2.2) becomes

(24.6.3.3)\[\frac{\Omega_P}{\Delta t} \sum_{i=1}^{N_s} v_{i,P}\left[(\rho^{*}Y_i)_P-(\rho^n Y_i^n)_P\right] + \frac{\Omega_P}{\Delta t}F_P \rho_P^{'} + \sum_{f\in\partial P} F_{f,P}\dot{m}_{f}^{*} + \sum_{f\in\partial P} F_{f,P}\dot{m}_{f}^{'} = b_{P,\mathrm{sp}}^{\mathrm{lag}}.\]

These are the characteristic normalization factors of the species-summed volumetric form. The closing code-mapping section lists the corresponding implementation stores.

24.6.4. Unsteady Contribution

Because \(Y_i^{'}\approx 0\), the implicit density-correction part is

(24.6.4.1)\[\sum_{i=1}^{N_s}\frac{\rho_P^{'}Y_{i,P}}{\rho_{i,P}} = F_P \rho_P^{'}.\]

Using the same EOS-compatible acoustic linearization as the mixture-based compressible forms,

(24.6.4.2)\[\rho_P^{'}\approx \frac{p_P^{'}}{a_P^2},\]

so the implicit unsteady term becomes

(24.6.4.3)\[\frac{\Omega_P}{\Delta t}F_P \rho_P^{'} \approx F_P \frac{\Omega_P}{a_P^2\Delta t} p_P^{'},\]

while the explicit lagged part is

(24.6.4.4)\[b_{P,\mathrm{time}}^{\mathrm{exp}} = -\frac{\Omega_P}{\Delta t} \sum_{i=1}^{N_s} v_{i,P} \left[(\rho^{*}Y_i)_P-(\rho^n Y_i^n)_P\right].\]

For BDF2, the same structure is retained: the usual time-integrator factors act on \(F_P\rho_P\Omega_P\) rather than only on \(\rho_P\Omega_P\).

24.6.5. Face Mass-Flux Correction

The corrected face mass flux is built in the same way as in the compressible mixture form. For a face between cells \(P\) and \(N\), oriented from \(P\) to \(N\),

(24.6.5.1)\[\dot{m}_f^{'} = \rho_f^{*}\dot{V}_f^{'} + \rho_f^{'}\dot{V}_{f,\mathrm{rel}}^{*},\]

where \(\dot{V}_{f,\mathrm{rel}}^{*}\) is the provisional relative volumetric flux.

Using the same pressure-difference correction and upwind acoustic linearization as in the mixture-based compressible form, write

(24.6.5.2)\[\dot{m}_f^{'} = -D_f\left(p_N^{'}-p_P^{'}\right) + K_{P,f}p_P^{'} - K_{N,f}p_N^{'}.\]

with

(24.6.5.3)\[D_f=\psi_f\,C_f, \qquad K_{P,f}= \chi_c \frac{\max\left(\dot{V}_{f,\mathrm{rel}}^{*},0\right)}{a_P^2}, \qquad K_{N,f}= \chi_c \frac{\max\left(-\dot{V}_{f,\mathrm{rel}}^{*},0\right)}{a_N^2}.\]

Equation (24.6.5.2) is just the split used in Pressure-Correction Equation on Unstructured Grids, specialized to the present form: \(D_f\) is the pressure-difference term and \(K_{P,f}\), \(K_{N,f}\) are the sign-dependent upwind acoustic additions.

Multiplying by the row-local species factor gives the contribution of one face to the \(P\) row:

(24.6.5.4)\[F_{f,P}\dot{m}_f^{'} = F_{f,P}\left(D_f+K_{P,f}\right)p_P^{'} - F_{f,P}\left(D_f+K_{N,f}\right)p_N^{'}.\]

This is the key algebraic step for the species-summed volumetric form. It shows directly how the usual pressure-correction machinery is being reused, but now every face contribution is weighted by a local species sum.

24.6.6. Final Coefficients

Collecting the terms in Eq. (24.6.3.3) gives

(24.6.6.1)\[a_P^{(\mathrm{sp})}p_P^{'} = \sum_{N\in\mathcal{N}(P)} a_{PN}^{(\mathrm{sp})}p_N^{'} + b_P^{(\mathrm{sp})}.\]

If the compressible additions are temporarily suppressed, the base neighbor coefficient from one face is

(24.6.6.2)\[a_{PN,f}^{(0)}=F_{f,P}D_f.\]

For one interior face between cells \(P\) and \(N\), the final facewise coefficients in the \(P\) row are therefore

(24.6.6.3)\[a_{P,f}^{(\mathrm{sp})}=F_{f,P}\left(D_f+K_{P,f}\right), \qquad a_{PN,f}^{(\mathrm{sp})}=F_{f,P}\left(D_f+K_{N,f}\right).\]

The total diagonal for row \(P\) is then

(24.6.6.4)\[a_P^{(\mathrm{sp})} = \sum_{f\in\partial P^{\mathrm{int}}} a_{P,f}^{(\mathrm{sp})} + \sum_{b\in\partial P^{p}} F_{b,P}D_b + F_P {}^{(p')}\!a_{P,\mathrm{time}},\]

with

(24.6.6.5)\[{}^{(p')}\!a_{P,\mathrm{time}} = \left(\beta_0-\beta_1\right) \chi_c \frac{\Omega_P}{a_P^2}.\]

Here \(\beta_0\) and \(\beta_1\) denote the current-state and lagged-state coefficients supplied by the selected time integrator.

The right-hand side is

(24.6.6.6)\[b_P^{(\mathrm{sp})} = -\sum_{f\in\partial P} F_{f,P}\dot{m}_{f}^{*} + b_{P,\mathrm{time}}^{\mathrm{exp}} + b_{P,\mathrm{sp}}^{\mathrm{lag}}.\]

This is the cleanest final coefficient statement for the species-summed volumetric branch: the species factors appear only because the underlying continuity-like equation has changed.

24.6.7. Post-Solve Corrections

Once Eq. (24.6.6.1) is solved, the correction step is the same as in Pressure-Correction Equation on Unstructured Grids:

(24.6.7.1)\[p_P = p_P^{*} + \alpha_p p_P^{'} ,\]

(24.6.7.2)\[\left(\mathbf{u}\right)_P = \left(\mathbf{u}^{*}\right)_P - \frac{\sum_f p_f^{'}\mathbf{A}_f}{{}^{(u)}\!a_P},\]

(24.6.7.3)\[\dot{m}_f = \dot{m}_f^{*} + \dot{m}_f^{'},\]

with \(\dot{m}_f^{'}\) given by Eq. (24.6.5.2). The \(F\) factors do not appear in these post-solve corrections; they appear only in the pressure-correction row because the governing continuity-like equation has changed.

For ideal-gas mixtures at common \((p,T)\),

(24.6.7.4)\[v_i=\frac{R_iT}{p} \;\Rightarrow\; F_P=\sum_i Y_i\frac{R_iT}{p} =\frac{T}{p}\sum_i Y_iR_i =\frac{1}{\rho},\]

so the species-based normalization collapses to the familiar mixture volumetric factor. In that limit, the species-summed volumetric form recovers the same normalization one would obtain by dividing mixture continuity by \(\rho\), but the implementation still reaches it through the explicit \(F\) factors.

24.7. Code Mapping Summary

The body of this appendix has used generic notation on purpose. The present section records how that notation maps back to the current implementation.

24.7.1. Family Names and Option Names

Mathematical family	Current option / path name	Practical meaning
constant-density mass form	default incompressible path, constantDensityPPrime	\(\sigma_p=1\), \(\sigma_\rho=1\)
face-density-divided form	incompressibleForm, variableDensityPPrime	\(\sigma_p=0\), \(\sigma_\rho=1\)
compressible bookkeeping path in incompressible mode	compressibleForm, constantDensityPPrime + compressiblePPrime bookkeeping	mass-form bookkeeping; acoustic additions remain inactive in incompressible mode
mixture mass form	default compressible path, compressiblePPrime	\(\sigma_p=1\), \(\sigma_\rho=1\) with acoustic additions
mixture volumetric form	compressibleVolumetricForm, compressiblePPrime	\(\sigma_p=1\), \(\sigma_\rho=0\)
species-summed volumetric form	compressibleSpeciesVolumetricForm	normalization carried by the factors \(F_P\), \(F_{f,P}\), and \(F_{f,N}\) rather than by setting \(\sigma_\rho=0\)

24.7.2. Symbol-to-Store Mapping

Analytic quantity	Current implementation name	Note
\(\sigma_p\)	pPrimeFactor	0/1 switch selecting whether the face density remains inside the pressure-difference coefficient.
\(\sigma_\rho\)	rhoNormFactor	0/1 switch selecting whether the assembled row is kept in mass form or divided by the cell density.
\(\chi_c\)	compressiblePPrimeFactor	Coefficient multiplying the implicit compressible density correction.
\(d_f\)	pPrimeCoefficient	Base face coefficient before face-density and rowwise normalization.
temporary row accumulators	pPrimeTemp_D, pPrimeTemp_B	Temporary diagonal and right-hand-side accumulators.
sign-absorbed stored row	pPrime_D, pPrime_L, pPrime_U	Stored coefficients after the sign convention of Eq. (24.1.2) has been applied.
\(F_P\)	F_C	Cellwise species-weighted normalization factor.
old-time factor associated with \(F_P\)	F_C_n	Used in the explicit unsteady term.
\(F_{f,P}\), \(F_{f,N}\)	F_L, F_R	Left and right facewise species-weighted normalization factors.
lagged species source contribution \(b_{P,\mathrm{sp}}^{\mathrm{lag}}\)	pPrimeTemp_B together with yProductionRate, yProductionRateSourceTerm, and ydiff	Chemistry and diffusion contributions enter through the lagged source assembly.
provisional relative volumetric flux \(\dot V_{f,\mathrm{rel}}^*\)	volFluxStar - gridVolFlux	Relative face flux used in the upwind acoustic correction.

The dedicated face-coupling rules for the species-summed volumetric form are currently constrained to SIMPLEC paths in the implementation, so in practice the factor \(\psi_f\) in Eq. (24.6.5.3) reduces to the SIMPLEC factor from Eq. (24.3.3.1) in that branch.