5. Momentum-Weighted Interpolation on Collocated Unstructured Grids#
The momentum equations provide cell-centered velocities. Continuity and convection require flow rates through the cell faces. Momentum-weighted interpolation uses the discrete momentum equation to construct these face fluxes while preserving their response to pressure differences between neighboring cells.
The stored area vector \(\vec A_f\) points from the left cell cl to the
right cell cr. The volumetric and mass fluxes use this same orientation:
Here \(\vec{u}_f\) is the face-averaged velocity. On a collocated grid, pressure and velocity are stored at cell centers. We estimate \(\vec{u}_f\) from the cell-centered momentum solution.
In a cell balance, the incidence sign \(\sigma_{Pf}\) is \(+1\) when \(P=\operatorname{cl}(f)\) and \(-1\) when \(P=\operatorname{cr}(f)\). An ordinary boundary face uses \(+1\) for its physical owner. These are the face-orientation conventions used below. On a fixed grid, continuity for cell \(P\) sums over its faces \(\mathcal F_P\):
The same mass flux carries momentum, energy, species, and other transported quantities through the face. For example, the convective flux of velocity component \(u_i\) is \(\dot m_f u_{i,f}\), where \(u_{i,f}\) is the reconstructed face value.
Momentum interpolation first constructs a provisional face flux \(\dot{\Omega}_f^*\) from the solved cell velocities. The star marks this provisional value. Pressure correction then adjusts the flux to satisfy continuity. The classical construction is also called Rhie–Chow interpolation [RC83].
The derivation starts with an orthogonal face and an unrelaxed momentum equation. Later sections add general face geometry, velocity under-relaxation, time history, balanced momentum sources, variable density, and moving grids. The face-flux formula also supplies the coefficient \(d_f\) used in the pressure-correction equation.
5.1. Available Formulations#
For the standard first- and second-order upwind fluxes,
momentumInterpolationOptions selects the formulation. The default setting,
momentumInterpolationOptions: <form=new>, selects Pascau interpolation.
The names are case-sensitive. The table links each formulation to its derivation.
Selection |
Formulation |
|---|---|
|
|
|
Pascau interpolation with consistent cell and face time history |
|
Classical interpolation with a Choi-style face-history correction |
|
Pascau refinement with a separately reconstructed acceleration source |
|
Refinement using the provisional velocity in that acceleration source |
|
Bartholomew–Denner interpolation with density and source weighting |
|
Alias of |
The conservativePressure modifier is described in
Consistency with Momentum Source Terms. Denner’s skewnessCorrection, highOrder, and
modified options are described in Geometry and Source Modifiers.
5.2. Why Cell-Velocity Interpolation Is Insufficient#
For a single interior face \(f\), label its cells in stored order: \(P=\operatorname{cl}(f)\) and \(N=\operatorname{cr}(f)\). The face area vector is
Here \(A_f\) is the face area and \(\hat{\vec{n}}_f\) is the stored unit normal. Interpolating the cell velocities with weights \(w_P\) and \(w_N\) gives
On a uniform one-dimensional mesh with spacing \(\Delta x\), let the cell pressures alternate between \(p_0+\Delta p\) and \(p_0-\Delta p\):
The cells on either side of cell \(i\) have the same pressure, so the centered pressure gradient is zero:
This pattern produces no cell pressure gradient and hence no pressure-gradient contribution to the cell velocity. The interpolated face flux cannot detect it. The pattern can persist while the flux satisfies continuity; this is called pressure checkerboarding.
The pressure difference between adjacent cells is
This difference captures the alternating pattern. Momentum-weighted interpolation uses it to make the face flux respond to that pattern. Starting from the discrete momentum equation, it removes the interpolated cell pressure-gradient term and replaces it with a term based on \(p_N-p_P\). The resulting face flux has the form
The next section derives this face-flux correction from the discrete momentum equation.
5.3. From the Momentum Row to the Face Flux#
Our goal is the face-flux formula outlined in Eq. (5.2.6): an interpolated cell-velocity flux plus a pressure-dependent correction. We first assemble the cell momentum equation, separate its pressure-gradient term, and then use those quantities to construct the face flux.
Start with a fixed grid, BDF1 time integration, and an unrelaxed momentum equation with no mass source. The velocity \(\vec u_P\) is the current cell velocity; the superscript \(n\) identifies the previous physical time level. In the cell-row assembly, \(P\) denotes the cell being assembled and \(\sigma_{Pf}\) supplies the sign of each incident face.
5.3.1. From the Transport Row to the Momentum Row#
Set \(\phi=u_i\) in the transport chapter’s algebraic row:
Here \({}^{(u_i)}\!a_P\) is the diagonal coefficient, \({}^{(u_i)}\!a_{P,N}\) multiplies the velocity in neighbor \(N\), and \({}^{(u_i)}\!b_P\) is the complete known right-hand side (RHS).
Use the componentwise viscous treatment developed in Chapter 4: the two-point diffusion term is implicit, and the remaining stress terms are explicit. With scalar implicit sources, all three velocity components share the same diagonal and neighbor coefficients. Denote these shared coefficients by \({}^{(u)}\!a_P\) and \({}^{(u)}\!a_{P,N}\), and collect the component right-hand sides into \({}^{(u)}\!\vec b_P\). The vector row is
To obtain the transport chapter’s continuity-added form, subtract \(\vec u_P\) times the discrete continuity equation from the momentum equation. The temporal part becomes
Expand the velocity difference to identify the temporal diagonal and known time-history term:
The coefficient of \(\vec u_P\) is \({}^{(u)}\!a_{P,\mathrm{time}}\). Moving the old-time term to the RHS gives the positive contribution \({}^{(u)}\!\vec b_{P,\mathrm{time}}\).
The same continuity subtraction changes each face’s convection contribution
to \(\sigma_{Pf}\dot m_f(\vec u_f-\vec u_P)\). Here \(\dot m_f\) retains
the stored cl-to-cr orientation. With first-order upwinding and \(N\)
denoting the cell across this face,
The factor \(\max(-\sigma_{Pf}\dot m_f,0)\) is the incoming mass-flux magnitude for \(P\). It equals \(\max(-\dot m_f,0)\) for \(P=\operatorname{cl}(f)\) and \(\max(\dot m_f,0)\) for \(P=\operatorname{cr}(f)\). Moving the neighbor term to the RHS adds this magnitude to both the diagonal and the neighbor coefficient. The implicit viscous flux \(\mathcal K_f(\vec u_N-\vec u_P)\) contributes \(\mathcal K_f\) to each coefficient in the same way. Here \(\mathcal K_f\) is the two-point diffusion coefficient in Eq. (4.2.1.24), evaluated with the momentum diffusivity.
Collecting the interior-face contributions gives
Here \(\mathcal F_P^{\mathrm{int}}\) is the set of interior faces of \(P\), and \(\mathcal F(P,N)\) contains the faces connecting \(P\) to \(N\). An integrated source \(-c_P\vec u_P+\vec s_P\) adds \(c_P\) to the diagonal and \(\vec s_P\) to the known RHS. Boundary conditions supply the corresponding boundary-face contributions.
The complete unrelaxed RHS is
The integrated pressure-gradient contribution is \(\vec{\mathcal P}_P\); its force in the momentum equation is \(-\vec{\mathcal P}_P\). Pressure is held at its current outer-iteration estimate during the momentum solve. The convection correction \({}^{(u)}\!\vec b_{P,\mathrm{conv}}\) contains the higher-order reconstruction terms. The viscous correction \({}^{(u)}\!\vec b_{P,\mathrm{diff}}\) contains the non-orthogonal and remaining stress terms. Both are evaluated from lagged reconstructed values. The term \({}^{(u)}\!\vec b_{P,\mathrm{src}}\) collects the known sources.
Contributions to the Momentum Row
The table collects the contributions to the momentum row in Eq. (5.3.1.2) for the unrelaxed BDF1 equation after continuity subtraction. Add the entries in each column to assemble the corresponding coefficient or RHS. The first-order convection and two-point viscous rows sum contributions from interior faces.
Contribution |
Adds to diagonal \({}^{(u)}\!a_P\) |
Adds to neighbor coefficient \({}^{(u)}\!a_{P,N}\) |
Adds to known RHS \({}^{(u)}\!\vec b_P\) |
|---|---|---|---|
BDF1 storage |
\(\rho_P^n\Omega_P/\Delta t\) |
\(0\) |
\((\rho_P^n\Omega_P/\Delta t)\vec u_P^n\) |
First-order convection |
\(\sum_{f\in\mathcal F_P^{\mathrm{int}}}\max(-\sigma_{Pf}\dot m_f,0)\) |
\(\sum_{f\in\mathcal F(P,N)}\max(-\sigma_{Pf}\dot m_f,0)\) |
\(0\) |
Two-point viscous flux |
\(\sum_{f\in\mathcal F_P^{\mathrm{int}}}\mathcal K_f\) |
\(\sum_{f\in\mathcal F(P,N)}\mathcal K_f\) |
\(0\) |
Higher-order convection correction |
\(0\) |
\(0\) |
\({}^{(u)}\!\vec b_{P,\mathrm{conv}}\) |
Explicit viscous correction |
\(0\) |
\(0\) |
\({}^{(u)}\!\vec b_{P,\mathrm{diff}}\) |
Integrated source \(-c_P\vec u_P+\vec s_P\) |
\(c_P\) |
\(0\) |
\(\vec s_P\) |
Pressure |
\(0\) |
\(0\) |
\(-\vec{\mathcal P}_P\) |
For the diagonal, each face sum covers all interior faces of \(P\). For the neighbor coefficient, it covers only the faces connecting \(P\) to that neighbor \(N\).
Boundary-face convection and viscous contributions are added according to the boundary condition. Prescribed boundary values enter the known RHS, while any coefficient multiplying the cell velocity enters the diagonal.
The diagonal and neighbor coefficients have units of mass per time. Every term in the vector row, including its RHS, has units of force.
These contributions supply the diagonal, neighbor coefficients, and known RHS in Eq. (5.3.1.2). We can now use that assembled equation to identify the pressure-gradient term in the cell velocity.
5.3.2. Collecting the Nonpressure Contribution#
To separate the pressure term, substitute the complete RHS from Eq. (5.3.1.7) into the vector momentum row. Collect the neighbor-velocity terms and the nonpressure RHS terms into \(\vec H_P\):
After solving the momentum equations, evaluate \(\vec H_P\) from the neighbor velocities and known RHS terms. The cell relation becomes
Fig. 5.3.2.1 The momentum row supplies a cell velocity. Momentum interpolation replaces its cell pressure-gradient term with a face pressure-difference term. Pressure correction then adjusts the face flux to satisfy continuity.#
For the classical gradient form,
Substitute this pressure contribution into the row and divide every term by \({}^{(u)}\!a_P\):
The second term is the pressure-gradient term in the cell velocity equation. Its coefficient, \(\Omega_P/{}^{(u)}\!a_P\), comes from dividing the momentum row by its diagonal. To see how this coefficient relates changes in the two quantities, hold \(\Omega_P\), \({}^{(u)}\!a_P\), and \(\vec H_P\) fixed. A change in the pressure gradient then gives
Here \(\Delta\) denotes a change. The pressure-gradient change is the input; the cell-velocity change is the output. This is the sensitivity of cell velocity to the pressure gradient for the fixed momentum row. The minus sign gives the direction of the velocity change, and \(\Omega_P/{}^{(u)}\!a_P\) sets its magnitude. We use this coefficient when constructing the pressure term at the face.
Define the momentum pseudovelocity as the nonpressure velocity contribution:
Rearrange Eq. (5.3.2.4) to evaluate it from the solved cell velocity:
We now have the pseudovelocity \(\hat{\vec u}_P\) and the coefficient \(\Omega_P/{}^{(u)}\!a_P\) multiplying the pressure gradient. The face construction uses both quantities from each adjacent cell.
5.3.3. Constructing the Flux on an Orthogonal Face#
To construct the face-normal velocity from the cell relation in Eq. (5.3.2.4), we need a face value of the pseudovelocity, a coefficient multiplying the pressure gradient, and an estimate of the face-normal pressure gradient.
For the face construction, use the stored ordering \(P=\operatorname{cl}(f)\) and \(N=\operatorname{cr}(f)\). Take an orthogonal face, whose cell-center line is parallel to its stored normal. First interpolate the nonpressure velocity from the two cells:
Interpolate the coefficient \(\Omega/{}^{(u)}\!a\) multiplying the pressure gradient using the same weights:
These interpolations supply the first two quantities. For the remaining quantity, the face-normal pressure gradient, use the pressure difference between the two cells. Let \(\vec x_P\) and \(\vec x_N\) be their center positions, with
The center spacing on this orthogonal face is
On this orthogonal face, \(\ell_f\) is also the center-to-center distance.
Approximate the face-normal pressure gradient by the pressure difference divided by that distance:
This gives the exact gradient for a linear pressure field. Equations (5.3.3.1), (5.3.3.2), and (5.3.3.5) now supply all three quantities. Using these three approximations in the face form of Eq. (5.3.2.4) gives
Multiply the normal velocity by the face area to obtain volume flow rate:
Substituting the reconstructed normal velocity defines the provisional flux:
5.3.4. Writing the Flux in Terms of Cell Velocities#
Equation (5.3.3.8) gives the provisional flux in terms of pseudovelocities. Substituting their definitions expresses the same flux using the solved cell velocities and identifies the momentum-interpolation correction. The two cell relations are
Substitute these expressions into the interpolated pseudovelocity term in Eq. (5.3.3.8) and collect the velocity and pressure terms:
The first term is the interpolated cell flux. The second adds back the cell pressure-gradient term, removing that term from the interpolated velocity. Substitute this expansion into Eq. (5.3.3.8):
Provisional Face Flux
This completes the construction outlined in Eq. (5.2.6). The first term interpolates the cell velocities; the remaining two terms form the momentum-interpolation correction.
Multiply the volume flux by the face density to obtain the mass flux:
5.3.5. A Compact Pressure-Difference Form#
Equation (5.3.4.3) contains two pressure terms: the positive term formed from the reconstructed cell gradients and the negative term formed from \(p_N-p_P\). We now combine them to express the correction as a difference between two estimates of the pressure change across the face.
To give these terms a common coefficient on the orthogonal face, interpolate the coefficient \(\Omega/{}^{(u)}\!a\) and the pressure gradient separately, then multiply the two interpolated quantities. This gives the product approximation:
Use \(\mathcal I_f(q)=w_Pq_P+w_Nq_N\) to write these interpolations more compactly. The center-to-center vector on this face is \(\vec r_{PN}=\ell_f\hat{\vec n}_f\). The component of the interpolated gradient normal to the face can therefore be written as
Apply the product approximation and substitute this geometric identity into the cell pressure term:
The two pressure terms now share a coefficient and denominator:
Factor out the common coefficient. The bracket compares the direct pressure change \(p_N-p_P\) with the change \(\mathcal I_f(\nabla p)\cdot\vec r_{PN}\) obtained from the cell gradients:
Substitute this combined pressure term into Eq. (5.3.4.3):
This is the compact form of the flux in Eq. (5.3.4.3), using the product approximation in Eq. (5.3.5.1). Its correction compares the adjacent-cell pressure difference with the pressure change obtained from the reconstructed cell gradients.
For a linear pressure field with exact reconstructed gradients, the two pressure changes agree and the correction is zero. For the alternating field in the checkerboard example, the cell gradients are zero, and the adjacent-cell pressure difference produces a nonzero face-flux correction.
The product approximation is the additional step that made this compact form possible. The difference between the original interpolated product and its approximation is
This approximation is exact when either the coefficient \(\Omega/{}^{(u)}\!a\) or the reconstructed pressure gradient is the same in both cells.
The classical Cartesian derivation develops the uniform-grid expression and its Taylor expansion, which describes the accuracy of the face interpolation for smooth pressure fields.
5.3.6. Geometry and the Two Pressure Reconstructions#
The preceding construction assumes that the cell-center line is parallel to the face normal. For a non-orthogonal face, we must account for their different directions when reconstructing the normal pressure gradient.
For a general face, use the transport chapter’s center-line geometry:
\(d_{PN}=\lVert\vec r_{PN}\rVert\),
\(\hat{\vec r}_{PN}=\vec r_{PN}/d_{PN}\), and
\(\ell_f=\vec r_{PN}\cdot\hat{\vec n}_f>0\), with \(P\) and \(N\) in the
stored cl-to-cr order. The angle between \(\hat{\vec r}_{PN}\) and \(\hat{\vec n}_f\)
measures non-orthogonality. Face-center skewness is the displacement of the
face centroid from the chosen interpolation point on the center line.
Fig. 5.3.6.1 Cell-center direction, face normal, and face-centroid displacement.#
Let \(\vec r_t=\vec r_{PN}-\ell_f\hat{\vec n}_f\) be the transverse part of the center-to-center vector. For a linear pressure field,
Divide by \(\ell_f\) and subtract the transverse term to obtain the normal pressure gradient. Using interpolated cell gradients for that term gives
The first term uses the adjacent-cell pressure difference. The second is the transverse correction, evaluated from lagged reconstructed cell gradients. Keeping only the first term gives the two-point reconstruction in Eq. (5.3.3.5) [FP02].
To use the corrected normal gradient in the expanded flux, Eq. (5.3.4.3), replace \((p_N-p_P)/\ell_f\) by the right-hand side of Eq. (5.3.6.2). This extends that flux construction to a non-orthogonal face. For a linear pressure field with exact reconstructed gradients, the cell-gradient and corrected face-gradient terms cancel.
The compact formula in Eq. (5.3.5.6) uses a second construction: retain the two-point face term \((p_N-p_P)/\ell_f\) and estimate the cell-gradient term along the center line:
On an orthogonal face, this is the exact identity in Eq. (5.3.5.2). On a non-orthogonal face, it is an additional approximation. Together with the product approximation in Eq. (5.3.5.1), this gives the compact formula, Eq. (5.3.5.6), on a non-orthogonal face. For a linear pressure field with exact reconstructed gradients, its two pressure changes still agree and the correction is zero.
The unrelaxed face-flux construction is now complete. The next section adds velocity under-relaxation to the momentum row and the corresponding face flux.
5.4. Classical Rhie–Chow Form and Under-Relaxation#
This section describes momentumInterpolationOptions: <form=old>, including
its Majumdar iteration-history term. The accepted form=CChoi alias selects
the same construction.
Iteration and Stage Notation
The superscript \(k\) identifies a completed outer, or Picard, iteration, while \(k+1\) identifies the iteration being advanced. A star identifies the provisional momentum-predictor stage before pressure correction, and a prime identifies the correction obtained from the pressure-correction equation. The basic algorithmic sequence is
A superscript \(n\) identifies a physical time level, independently of \(k\). When only the predictor stage matters, \(\vec u^*\) and \(\dot\Omega_f^*\) stand for \(\vec u^{*,k+1}\) and \(\dot\Omega_f^{*,k+1}\).
5.4.1. Relaxation of the Momentum Row#
Let \({}^{(u)}\!a_P\) and \({}^{(u)}\!\vec b_P\) denote the unrelaxed assembled row, and let \(0<\alpha_v\leq1\) be the velocity relaxation factor. Add \([(1-\alpha_v)/\alpha_v]{}^{(u)}\!a_P (\vec u_P^*-\vec u_P^k)\) to the left side of that row. This addition vanishes when the iteration stops changing. Moving its known \(\vec u_P^k\) term to the right gives
The diagonal is divided by \(\alpha_v\); the neighbor coefficients are unchanged. The last term is a known source that restores the original row when the iteration converges. In the algorithm chapter’s notation,
Dividing Eq. (5.4.1.1) by its diagonal, and collecting the nonpressure neighbor and source terms into \(\vec H_P\), gives the velocity relation below. Neighbor velocities inside \(\vec H_P\) are those of the implicit predictor solution, not a separately relaxed explicit update.
The resulting provisional cell velocity is
The bracket is the same momentum relation as Eq. (5.3.2.4), evaluated with the new predictor’s neighbor velocities. With \(\alpha_v=1\), the previous-iteration term disappears. With \(\alpha_v<1\), part of the previous corrected velocity is retained.
The volumetric flux obtained by interpolating the previous cell velocities is
This interpolation generally differs from the stored face flux \(\dot{\Omega}_f^{\,k}\). Their difference is
The difference is already known from iteration \(k\). It contains the face-specific contribution that is not reproduced by simply interpolating the cell velocities. Majumdar’s construction carries the fraction \(1-\alpha_v\) of this previous face-flux difference into the new provisional flux [Maj88].
The new predictor already contains the interpolated part of the previous
cell velocity. Carry the remaining previous face-flux difference separately.
Combining it with the two pressure terms gives the expanded form=old result:
Each coefficient \(\Omega/{}^{(u)}\!a\) uses the complete unrelaxed diagonal. Integrated pressure and body-source replacements are developed in Consistency with Momentum Source Terms.
At convergence, the predictor and corrected quantities agree between iterations. Denote the remaining difference between the face flux and the interpolated cell flux by
Subtract the interpolated velocity term from Eq. (5.4.1.6). The same difference now occurs on both sides:
Move the last term to the left. Its coefficient becomes \(1-(1-\alpha_v)=\alpha_v\). Dividing by \(\alpha_v>0\) leaves
This is the unrelaxed pressure contribution. The previous face-flux difference is therefore what removes the explicit relaxation factor from the converged relation. Its temporal dependence is addressed next.
We now combine the two pressure terms in the expanded flux, Eq. (5.4.1.6), using the steps from A Compact Pressure-Difference Form. Interpolate the coefficient and pressure gradient separately, then multiply them. Use the interpolated gradient to estimate the pressure change between cell centers. These substitutions give the compact flux formula:
The history term is unchanged. The pressure terms carry the additional approximations used to obtain the compact form.
Let \(\vec x_f\) be the face-center position. For this family, the interpolation weights are based on the cell-center-to-face distances:
5.4.2. Relaxing Only the Spatial Diagonal#
An alternative is to relax only the non-temporal part of the row. Write \(\widetilde a_P={}^{(u)}\!a_{P,\mathrm{spatial}}\) and let \(T_P={}^{(u)}\!a_{P,\mathrm{time}}\). Whole-diagonal relaxation and spatial-only relaxation then give, respectively,
The deferred relaxation source is respectively
\((1-\alpha_v)(\widetilde a_P+T_P)\vec u_P^k/\alpha_v\) or
\((1-\alpha_v)\widetilde a_P\vec u_P^k/\alpha_v\).
Their different predictor diagonals give different coefficients multiplying
the pressure gradient in the velocity equation at finite iteration. The
derivations here use whole-diagonal relaxation. Spatial-only relaxation is
an algebraic alternative, not another form selector in the table above.
5.5. Time-Consistent Momentum Interpolation#
The Pascau construction in this section corresponds to
momentumInterpolationOptions: <form=new>. The complete one-step derivation
uses a fixed grid and the continuity-subtracted BDF1 momentum row. Time-step
independence below refers to the converged steady limit; transient solutions
still have the truncation error of their time discretization.
The classical momentum-interpolation construction uses the assembled momentum diagonal in the coefficient multiplying the cell pressure gradient. In a transient calculation, that diagonal contains a contribution from the time derivative as well as contributions from convection, diffusion, and any other terms treated implicitly. The temporal contribution can be seen directly from the discretized momentum storage term.
For a representative BDF1 step, after the continuity-consistent velocity form has been obtained and with \(\rho_P\Omega_P\) evaluated at the previous physical time level, the temporal contribution is
The current physical-time index on \(\vec{u}_P^{*,k+1}\) is suppressed; the star and \(k+1\) identify the momentum predictor within that time step. When the momentum row is assembled, the coefficient multiplying \(\vec{u}_P^{*,k+1}\) is placed on the left-hand side. The known term containing \(\vec{u}_P^n\) is moved to the right-hand side. Thus BDF1 adds \(\rho_P^n\Omega_P^n/\Delta t\) to the diagonal of every velocity-component row.
Within this continuity-subtracted temporal split, let \(\lambda_0\) denote the coefficient of the current-time contribution retained implicitly. BDF1 has \(\lambda_0=1/\Delta t\); the multistep split is developed below. The complete momentum diagonal is then
The shorthand \(\widetilde a_P={}^{(u)}\!a_{P,\mathrm{spatial}}\) is what remains after the current-time temporal contribution is removed. It contains the convective, diffusive, and other non-temporal contributions treated implicitly in the momentum row.
The temporal coefficient in this split belongs to the stated momentum-row form. It is not a universal replacement for the conservative or multistep temporal diagonal. Also, \(\widetilde a_N\) is the non-temporal diagonal of cell \(N\), not the coupling \({}^{(u)}\!a_{P,N}\) in the row for \(P\).
Using the complete \({}^{(u)}\!a_P\) in the classical coefficient \(\Omega_P/{}^{(u)}\!a_P\) introduces a temporal contribution into the face interpolation. The coefficient depends on the time step:
For BDF1, reducing \(\Delta t\) increases \(\lambda_0=1/\Delta t\) and therefore reduces \(\Omega_P/{}^{(u)}\!a_P\) when the other quantities are held fixed. If no matching face-flux history is carried, a converged nominally steady result can retain a residual dependence on the time step used to reach it. Time-consistent formulations use \(\Omega_P/\widetilde a_P\) and carry the corresponding physical-time history at the same face location as the new flux [Cho99, Pas11].
The ratio of the temporal diagonal to the non-temporal diagonal is
For BDF1, this becomes \(\delta_P=(\rho_P^n\Omega_P^n/\Delta t)/\widetilde a_P\). The ratio therefore measures the temporal contribution to the diagonal relative to the non-temporal contribution. The complete diagonal can be written
For the face shared by cells \(P\) and \(N\), the same geometric weights used in the momentum interpolation give
5.5.1. Deriving the Complete One-Step Face Relation#
The nonpressure group \(\vec H_P\) contains the old-time momentum as well as the spatial terms. Expose those two pieces before dividing the row:
Divide the spatial group by \(\widetilde a_P\) and call the resulting velocity \(\vec h_P\):
Thus \(\vec H_P=\widetilde a_P\vec h_P+ {}^{(u)}\!a_{P,\mathrm{time}}\vec u_P^n\). Substitute this and \({}^{(u)}\!a_P=\widetilde a_P+{}^{(u)}\!a_{P,\mathrm{time}}\) into the relaxed momentum row:
Multiply by \(\alpha_v/\widetilde a_P\). The two ratios that appear are
Replacing these ratios by their names gives
For the face, interpolate \(\vec h\), \(\delta\), and \(\Omega/\widetilde a\) with the selected geometric weights. Use the face’s own old-time and old-iteration fluxes in its history terms. The analogous face relation is defined by
This is an interpolation model for the face, not a new equation solved on a face control volume. Eliminate \(\mathcal I_f(\vec h)\) using the two cell rows. Rearranging and interpolating Eq. (5.5.1.5) gives
The ratios remain inside the interpolation of each cell velocity. For example, the provisional-velocity term expands as
It reduces to ordinary interpolation if \(\delta_P=\delta_N\). For unequal ratios, each cell velocity retains its own multiplier.
Now multiply Eq. (5.5.1.7) by \(A_f\), project onto the face normal, and substitute it into the face relation. Collect terms from the same time or iteration level:
The positive history terms came from the face relation. The negative ones came from eliminating the cell quantity \(\vec h\). Their subtraction replaces the history already carried by the cell velocities with the face’s own history.
5.5.2. Face-Flux History Terms#
Call the last bracket \(\eta_f^n\), the weighted difference between the previous-time face flux and the flux represented by the cell velocities:
Call the previous-iteration bracket \(\eta_f^k\):
After division by \(1+\delta_f\), the two history groups contribute
The cell and face histories are therefore treated as one discrete temporal system. With these definitions, the full one-step result is
Complete Pascau Construction
The four terms are the weighted provisional velocity, removal of the cell pressure-gradient term, insertion of the face pressure-difference term, and matched iteration/time history. Changing only the pressure coefficient does not produce this construction. A conservative pressure replacement changes the cell pressure term as described in Consistency with Momentum Source Terms.
5.5.3. Reading the Pressure Coefficient from the Flux#
In the completed face formula, hold the predictor, histories, gradients, and coefficients fixed. Only the direct pressure-difference term changes when \(p_N-p_P\) changes. Multiply its volume-flux coefficient by \(\rho_f\) to obtain the mass-flux correction coefficient \(d_f\). The signed sensitivity of face mass flux to \(p_N-p_P\) is \(-d_f\); increasing that pressure difference decreases the flux in the stored direction.
The classical old selection uses the complete momentum diagonal in
its mass-flux correction coefficient, while the one-step time-consistent
new selection uses the non-temporal diagonal together with
\(1+\delta_f\):
The superscripts \(\mathrm{old}\) and \(\mathrm{new}\) identify the Stream
momentum-interpolation selections; they are not physical-time levels. Both
expressions give the magnitude of the face mass-flux sensitivity to the
adjacent-cell pressure difference, with the remaining predictor quantities
held fixed. The new coefficient makes two changes: it uses the non-temporal
cell ratios \(\Omega/\widetilde a\), and it divides the resulting face
coefficient by \(1+\delta_f\).
5.5.4. Checking the Steady Limit#
To check the steady limit, define the spatial pressure contribution
At a converged steady state, the old-time, old-iteration, and provisional velocities all equal \(\vec u\), and the corresponding face fluxes equal \(\dot\Omega_f\). Insert these equalities into the expanded substitution. Move the history terms containing face flux to the left, and collect the cell-velocity terms on the right:
Both brackets reduce to \(\alpha_v\), even when the two cell ratios differ. Dividing by \(\alpha_v>0\) gives
Neither the temporal ratio nor the relaxation factor remains in this steady-limit relation. This cancellation uses the weighted predictor as well as both history terms; the history terms alone do not vanish separately for unequal cell ratios.
5.5.5. Multistep and Crank–Nicolson Extensions#
For constant step and fixed cell mass \(M_P=\rho_P\Omega_P\), start with the usual BDF2 storage term. Regroup its three velocity values as two increments:
The first current-velocity term is implicit; the second current-velocity term is lagged during assembly. At convergence this is the usual BDF2 derivative. It gives \(\delta_P=2M_P/(\Delta t\widetilde a_P)\) and \(\epsilon_P=M_P/(2\Delta t\widetilde a_P)\). The fully implicit BDF2 diagonal \(3M_P/(2\Delta t)\) must not be substituted for the leading split diagonal while retaining these history terms.
In the face construction, the older increment contributes a current-velocity term as well as an old-time term. Collecting the current terms gives the face multiplier
Here \(\epsilon_f=\mathcal I_f(\epsilon)\) uses the older-history coefficient just defined; it is not a stored value of \(\delta\).
With \(\Xi_P=1+\delta_P-\alpha_v\epsilon_P\) and \(\Xi_f=\mathcal I_f(\Xi)\), the corresponding assembled face result is
The definitions of \(\eta_f^k\) and \(\eta_f^n\) retain the current \(\delta\) coefficients. The face construction evaluates the split current-level contribution with the predictor velocity; the cell momentum assembly may lag that contribution during its solve. This additional face choice agrees with the cell split at outer-iteration convergence. This formula displays all three changes from BDF1: the predictor weighting, denominator, and older-history difference. For variable steps, with \(h=\Delta t_n\) and \(k_t=\Delta t_{n-1}\), the corresponding split coefficients are
Each history carries its own mass. Subtracting the two current-level coefficients recovers the usual leading BDF2 coefficient at constant mass. The variable-step derivative and continuity subtraction are developed in (10.3.2.6) and (10.3.3.1). The time-step symbol \(k_t\) here is unrelated to the outer-iteration index \(k\).
Crank–Nicolson can be normalized by dividing its half-weighted current spatial residual by \(1/2\). The temporal coefficient is then \(2M_P/\Delta t\) and the known old spatial residual also enters the RHS. This has the one-step face-history structure above with its correctly assembled \(\widetilde a_P\) and \(\delta_P\); it has no \(n-1\) history. In particular, \(\Delta t\delta_P/\rho_P=2\Omega_P/\widetilde a_P\) in this normalization. Replacing \(\Omega/\widetilde a\) by \(\Delta t\delta/\rho\) without that factor of two would change the coefficients multiplying the pressure terms. See (10.3.2.11) for the temporal normalization.
5.5.6. Choi-Style Face History#
The selection momentumInterpolationOptions: <form=Choi> retains the
classical relaxed flux and adds the face-history term below. It is a
Choi-style construction with a common face ratio in the history correction;
the original literature contains different cell/face weighting choices.
Choi’s unsteady correction retains the face’s own previous-time flux rather than recreating its history solely by interpolating previous cell velocities. In the one-step form, the characteristic correction is
Pascau’s construction retains the cell-specific \(\delta_P\) and \(\delta_N\) inside the interpolation, as in Eqs. (5.5.2.1) and (5.5.2.2), and provides a systematic criterion for determining whether a face-velocity formula remains dependent on relaxation or time-step size [Pas11].
The main distinctions within this family are therefore:
the classical form carries nonlinear-iteration history but not a consistently derived physical-time face history;
the Choi form carries an explicit face-flux history;
the Pascau form separates the temporal diagonal and constructs the face and cell histories with corresponding \(\delta\) factors; and
later Pascau-style refinements reorganize the time-history difference so that it remains compatible with balanced acceleration and body-force states.
These variants change the construction of the provisional flux and its mass-flux correction coefficient. They do not change the continuity equation that consumes the flux.
For form=Choi, the complete one-step pressure/velocity formula is Eq.
(5.4.1.6) plus \(\mathcal C_{f,\mathrm{Choi}}\). Its mass-flux correction
coefficient still uses the complete momentum diagonal. For form=new, use
Eq. (5.5.2.4), which also changes the predictor and
iteration-history weighting. Equal temporal ratios make several factors
coincide; unequal ratios reveal why these two constructions are different.
5.6. Consistency with Momentum Source Terms#
The conservativePressure modifier uses the integrated cell pressure
contribution in the old, new, Choi, and newer formulations. For example,
momentumInterpolationOptions: <form=new, conservativePressure> enables this
replacement for Pascau interpolation. It changes the pressure term removed
from the cell velocity and retains the direct adjacent-cell pressure difference.
Denner uses the integrated contribution as part of its formulation. The
newest formulation uses the reconstructed cell pressure gradient.
Equation (5.3.2.3) represents the integrated pressure-gradient contribution as \(\Omega_P\nabla p_P\). The actual momentum assembly may instead evaluate that contribution from face pressures,
The incidence sign gives each stored face contribution its sign in the cell row. The reconstruction of \(p_f\) determines the discrete pressure force. The row may also contain body-force, non-inertial, or surface-force contributions that balance pressure in an equilibrium state. If momentum interpolation removes \(\Omega_P\nabla p_P\) from the cell velocity while the momentum equation actually applied \(\vec{\mathcal{P}}_P\), the face construction no longer represents the discrete momentum equation that produced the cell velocity.
The consistent construction uses the actual integrated source:
and interpolates \(\vec{\mathcal{P}}_P/{}^{(u)}\!a_P\) when removing the pressure contribution already contained in the cell velocities. The same principle applies to a discrete body-force source \(\vec{\mathcal{B}}_P\).
For a balanced state,
the face-pressure and body-force corrections should cancel when their direct face reconstructions also satisfy the corresponding balance. This requires compatible interpolation and geometry at both locations. It is essential for hydrostatic, free-fall, rotating-frame, and other source-balanced solutions. A mismatch can generate a nonzero face flux even when the cell momentum residual is identically zero.
Discrete Compatibility Condition
Preserving an equilibrium requires both cancellation of the integrated cell forces and compatible direct pressure/source reconstructions at the face. Cell cancellation alone is insufficient. These are discrete balance conditions, beyond agreement in the continuum limit.
5.6.1. Removing and Replacing a Source Contribution#
Let \(\vec S_P\) be a body force per unit volume, and let \(\vec{\mathcal B}_P\) be its integrated contribution to the cell momentum row. In the divided Pascau row, this source contributes
Therefore the weighted predictor in the face formula already carries \(\alpha_v A_f\mathcal I_f(\vec{\mathcal B}/\widetilde a) \cdot\hat{\vec n}_f\). A separate face source must replace this contribution.
To construct that face source, start with the pressure change needed to balance it. Integrate along the two segments from \(P\) to the face center and from the face center to \(N\). With a constant source on each segment,
The line estimate \(L_f(\vec S)\) has pressure units. Dividing by \(\ell_f\) gives the direct face source used in the normal-velocity relation. The Pascau flux therefore needs the following change:
Call this pair \(\mathcal C_{S,f}(1+\delta_f)\). Other formulations use the same removal and replacement with a different face denominator. Writing that denominator as \(\Xi_f\) gives the reusable expression
For Pascau, \(\Xi_f=1+\delta_f\). For the classical formulation, replace \(\widetilde a\) by the complete diagonal and use \(\Xi_f=1\).
The pressure contribution with integrated cell pressure is
Add the pressure and source corrections and collect their cell and face terms separately:
The first line vanishes when the integrated forces balance in each cell. The second vanishes when the adjacent pressure difference equals the source line estimate. Both are required to preserve the equilibrium. Implicit source contributions remain in the momentum diagonal; include only the corresponding explicit remainder in the source replacement. The porous-media appendix develops that split.
5.6.2. A Hydrostatic Face#
For constant gravity, let \(\vec S=\rho\vec g\) and take an orthogonal face. Let the distances from the two centers to the face be \(h_P\) and \(h_N\). If density is piecewise constant on the two segments, hydrostatic balance gives
For example, \(h_P=h_N=0.5\ \mathrm m\), \(\rho_P=1\ \mathrm{kg/m^3}\), \(\rho_N=1000\ \mathrm{kg/m^3}\), and a downward unit normal with \(\vec g\cdot\hat{\vec n}_f=10\ \mathrm{m/s^2}\) give a pressure rise of \(5005\ \mathrm{Pa}\). If the cell pressure and gravity forces are also assembled to cancel, the two face corrections cancel exactly and a stationary predictor produces zero flux. A harmonic density multiplied by the full spacing would give a different pressure rise here; the density average must follow the particular term being reconstructed.
For a static interface with surface tension, the same requirement is \(\vec{\mathcal P}=\vec{\mathcal B}_{\sigma}\) at the cells and a matching pressure jump at the face. Curvature and interface reconstruction determine whether this balance is achieved. The interpolation principle alone does not guarantee that an approximate curvature field is in equilibrium.
5.7. Separately Reconstructed Acceleration#
The selections newer and newest change the treatment of the unsteady
contribution in new. The following derivations use BDF1 on a fixed grid;
they do not prescribe a multistep extension. Let
These have force-per-volume units. They enter as negative sources when the storage term is moved from the left to the right of the momentum balance. Their integrated cell contribution is \(\Omega_P\vec T_P\). Use Eq. (5.6.1.4) to define \(\mathcal C_{T,f}\) with \(\vec S=\vec T\) and \(\vec{\mathcal B}=\Omega\vec T\).
5.7.1. Acceleration from the Previous Iteration#
This is momentumInterpolationOptions: <form=newer>. Rewrite the temporal
and iteration terms in Eq. (5.5.1.5) as
This step adds and subtracts \(\alpha_v\delta_P\vec u_P^k\). The last difference is the previous-iteration acceleration because \(\delta_P=\rho_P^n\Omega_P/(\Delta t\widetilde a_P)\). The cell row can now be solved for its spatial nonpressure contribution:
At the face, use the previous face flux for the iteration term and the line estimate \(L_f(\vec T^k)/\ell_f\) for the acceleration. Its sign is negative because storage has been moved to the source side of the row:
The last line applies the body-source replacement already derived. Substitute the two cell expressions for \(\alpha_v\vec h\). The pressure terms reproduce \(\alpha_v\mathcal C_{p,f}\). The acceleration terms are
The remaining predictor and iteration terms carry straight through. Collecting them with these pressure and source pairs gives
Here \(\mathcal C_{p,f}\) is the spatial pressure expression in Eq.
(5.5.4.1). With conservativePressure, replace
\(\alpha_v\mathcal C_{p,f}/(1+\delta_f)\) by
\(\mathcal C_{p,f}^{\mathrm{int}}(1+\delta_f)\) in the divided flux equation.
The pressure denominator remains \(1+\delta_f\). The changed terms are the
history combination and the explicit acceleration reconstruction.
5.7.2. Acceleration from the Provisional Velocity#
This is momentumInterpolationOptions: <form=newest>. The predictor
acceleration satisfies
Insert this identity into the Pascau cell row. Move \(\alpha_v\delta_P\vec u_P^*\) from the right to the left. The coefficient of \(\vec u_P^*\) becomes \((1+\delta_P)-\alpha_v\delta_P=1+(1-\alpha_v)\delta_P\):
Call this coefficient \(\Xi_P=1+(1-\alpha_v)\delta_P\), and interpolate it to obtain \(\Xi_f=1+(1-\alpha_v)\delta_f\). Rearranging the cell row gives
Construct the face relation using the same coefficient and the predictor acceleration line estimate:
After substituting the cell expression for \(\alpha_v\vec h\), the terms group as follows:
Replace the two underbraced groups by their names to obtain
The predictor weighting and the pressure and source terms all use this new denominator. The direct pressure coefficient is therefore
For this coefficient, hold the computed predictor and its acceleration fixed while changing the adjacent-cell pressure difference.
In uniform free fall, \(\vec T=\rho\vec g\) when the velocity increment
matches the time-discrete acceleration. A matching body-source and
acceleration reconstruction then cancels their face corrections. For
newer, the previous iterate must already represent that increment; for
newest, the acceleration uses the new predictor. This explains the intended
balance mechanism. Convergence and preservation of a particular discrete
state also depend on the momentum solve, source assembly, and chosen
iteration, so this identity is not a blanket robustness guarantee.
5.8. Denner Unified Momentum-Weighted Interpolation#
This section describes momentumInterpolationOptions: <form=Denner>.
The Bartholomew–Denner construction uses the same basic operation as the
preceding sections: remove the pressure and source contributions carried by
the cell velocities and insert their direct face counterparts. Its distinct
choices concern the shared face coefficient, density weighting, and face
history [BDAA+18].
5.8.2. Removing Cell Pressure and Sources, Then Inserting Face Values#
For pressure, the direct face value is \((p_N-p_P)/d_{PN}\). The cell value is obtained from the integrated pressure forces, divided by cell mass, interpolated, and multiplied by \(\rho_{h,f}^k\). Their difference is
Thus \(-\mathcal D_f A_f\mathscr R_{p,f}\) adds back the pressure term represented by the cells and subtracts the face pressure-difference term, with the same signs as in the classical construction.
For a body source, the cell contribution enters velocity with a positive sign. Its replacement therefore has the opposite order: remove the cell contribution and add the direct face source. Define
Both differences have pressure-gradient units. Their combined spatial contribution to volume flux is \(\mathcal D_f A_f(-\mathscr R_{p,f}+\mathscr R_{b,f})\). The term residual here means a difference between two reconstructions; it is not the residual of a separate equation solved on the face.
For the gravity source, the direct face estimate uses the complementary weights
Their role follows from the source line estimate. On an orthogonal face, let \(h_P\) and \(h_N\) be the distances from each center to the face. Then \(w_N=h_P/d_{PN}\) and \(w_P=h_N/d_{PN}\), giving
This is the pressure rise obtained in the hydrostatic example. The velocity weights would exchange the two segment lengths in this source estimate. For gravity with a consistently subtracted reference density, \(\vec S_P=(\rho_P-\rho_{\mathrm{ref}})\vec g\).
The integrated cell forces must use the matching reconstruction. With \(p_f=w_Pp_P+w_Np_N\), the pressure-gradient contribution is \(\vec{\mathcal P}_P=\sum_{f\in\mathcal F_P}\sigma_{Pf}p_f\vec A_f\). The corresponding interior-face additions to the body source are
Both expressions use the same stored area vector. Reversing the center line and area together leaves the physical source unchanged. At a physical boundary, the one-sided addition is \([\vec S_P\cdot(\vec x_f-\vec x_P)]\vec A_f\) for the physical owner. These source contributions are the ones removed in Eq. (5.8.2.2).
If \(\vec{\mathcal P}=\vec{\mathcal B}\) in the cells and \((p_N-p_P)/d_{PN}=\vec S_f\cdot\hat{\vec r}_{PN}\) at the face, the two residuals are equal. The negative pressure correction and positive source correction then cancel. This checks both the signs and the density weighting.
5.8.3. Including Face Time History#
First take \(\alpha_v=1\). As in the classical derivation, separate the face flux from the flux represented by the cell velocities. Write the current difference as \(E_f\) and the stored differences as \(E_f^j\):
Using the shared face coefficient, the face time coefficient is \(c_f=\rho_{h,f}^n\lambda_0\). Subtracting the interpolated cell storage increment from the face storage increment gives
Storage is on the left of the momentum balance, so its contribution on the source side is \(-c_f(E_f-E_f^n)\). Combine it with the pressure and body-source differences. Let \(\beta_f\) be the common geometric multiplier on those two spatial terms; the default is \(\beta_f=1\), and the optional angular corrections are given below. The resulting face model is
The choice of a shared face coefficient and harmonic density enters in forming this model. The next steps are algebra. Move the term containing the unknown current difference to the left:
Divide by \(1+c_f\mathcal D_f\) and label the two resulting factors:
Thus the coefficient multiplying the pressure and source terms is \(F_f\), and the old-time weight is \(G_f\):
For under-relaxation, retain the fraction \(1-\alpha_v\) of the previous iteration’s flux difference and use \(\alpha_v\) of the new difference just derived. This is the same face-history construction used for the classical form:
Add back the interpolated predictor and expand the two stored differences:
For equal density, equal cell ratios \(\Omega/\widetilde a\), midpoint weights, and an orthogonal face, \(c_f\mathcal D_f=\delta_f\). The Denner and Pascau history and pressure factors then coincide. Equal spacing alone is insufficient when density or the ratio \(\Omega/\widetilde a\) differs between the cells.
At a converged steady state, all three face–cell flux differences coincide. Subtract the iteration-history term and divide by \(\alpha_v\) to get \((1-G_f)E_f=\beta_f F_f A_f(-\mathscr R_{p,f}+\mathscr R_{b,f})\). Since \(1-G_f=1/(1+c_f\mathcal D_f)\), the ratio \(F_f/(1-G_f)=\mathcal D_f\). The complete steady flux difference therefore uses the spatial coefficient \(\mathcal D_f\), with no remaining temporal factor.
5.8.4. Geometry and Source Modifiers#
The factor \(\beta_f\) scales the pressure and source terms in the face flux
when the center line and normal are not parallel. The setting skewnessCorrection=1 uses \(\beta_f=1\);
this is also the default. The choices skewnessCorrection=2 and
skewnessCorrection=3 use, respectively,
These factors depend on the angle between the face normal and the center line. All three coincide on an orthogonal face. The third changes \(A_f/d_{PN}\) into \(A_f/\ell_f\); the same factor multiplies both pressure and source residuals.
With highOrder, define
\(\vec r_f=\vec x_f-(w_P\vec x_P+w_N\vec x_N)\) and reconstruct the
pressure used in the integrated cell force as
This alters the cell pressure contribution and hence the explicit residual. When these reconstructed cell gradients are held fixed during pressure correction, the direct adjacent-cell coefficient and matrix connections are unchanged.
With modified, the direct gravity/source estimate instead uses
\(L_f(\vec S)/d_{PN}\) in place of \(\vec S_f\cdot\hat{\vec r}_{PN}\).
The integrated cell source additions in Eq.
(5.8.2.5) replace
\(\vec S_f\cdot\vec r_{PN}\) by \(L_f(\vec S)\) at the same time.
The pressure and cell-source removal remain projected along the center line.
Both the direct source and the cell source must be changed together.
5.8.5. Interior and Boundary Pressure Coefficients#
Only the direct face-pressure difference contains the unknown neighboring pressure correction. Linearizing that term gives the Denner interior-face coefficient
At a boundary face there is one adjacent cell. With
the corresponding coefficient is
For a specified-pressure boundary, this coefficient contributes to the pressure-correction diagonal. For a boundary whose mass flux is prescribed, the pressure correction does not alter that flux, so no corresponding unknown face-pressure correction is assembled.
5.9. Connection to the Pressure-Correction Equation#
The face-flux formula gives a second sensitivity relation: the input is the pressure difference between adjacent cells, and the output is face mass flux. To measure this sensitivity, vary the pressure difference while holding the computed predictor, history, cell gradients, sources, density, and momentum coefficients fixed. Define the mass-flux correction coefficient \(d_f\) as the negative of this signed sensitivity:
With \(P=\operatorname{cl}(f)\) and \(N=\operatorname{cr}(f)\), a positive \(d_f\) means that increasing \(p_N-p_P\) reduces the flux in the stored direction. Thus \(d_f\) gives the magnitude of the mass-flux change per unit change in the pressure difference. The corresponding volume-flux correction coefficient is \(d_f/\rho_f\). These coefficients retain the local relation between pressure and flux for use in the continuity equation.
All of the preceding interpolation families supply the same base linear pressure-difference correction. For an interior face with \(P=\operatorname{cl}(f)\) and \(N=\operatorname{cr}(f)\), the volumetric-flux and mass-flux corrections for iteration \(k+1\) with SIMPLE are
The pressure-correction difference is the input to these relations; the resulting flux corrections are the outputs. Substituting the mass-flux relation into continuity gives the pressure-correction equation. Solving that equation determines the pressure corrections that satisfy linearized continuity, using the coefficient \(d_f\) from the face-flux formula.
Adding these corrections to the provisional fluxes defines the corrected fluxes that leave the outer iteration:
These equalities describe the algorithmic correction step; they are not interpolation approximations. The corrected flux is carried into the next outer iteration’s convection assembly. Compressible density-flux corrections are additional terms developed in the pressure-correction chapter; they are not included in the fixed-density flux correction above.
5.9.1. One Face in the Pressure Row#
For the row of an arbitrary cell \(P\), let \(N\) be its neighbor across face \(f\). The face contribution to fixed-density steady continuity is \(\sigma_{Pf}\dot m_f^*+d_f(p_P'-p_N')\). Collecting the pressure unknowns into \({}^{(p')}\!a_Pp_P'=\sum_N{}^{(p')}\!a_{P,N}p_N'+{}^{(p')}\!b_P\) gives
The neighbor has \(\sigma_{Nf}=-\sigma_{Pf}\). Its diagonal and neighbor contributions are again \(d_f\), and its RHS contribution is \(+\sigma_{Pf}\dot m_f^*\). The two cell contributions conserve the face mass flow. Temporal continuity defects, density corrections, and normalization add their own terms as derived in the next chapter.
At a specified-pressure boundary, a known boundary correction \(p_{\mathrm{bnd}}'\) gives \(d_{\mathrm{bnd}}p_P'\) on the left and \(d_{\mathrm{bnd}}p_{\mathrm{bnd}}'\) on the right. If the predictor already satisfies the prescribed pressure, \(p_{\mathrm{bnd}}'=0\). A prescribed-flux boundary contributes its known flux to the continuity defect and has no unknown pressure-induced flux term. For the classical/Pascau boundary flux correction, the one-sided coefficient uses the adjacent cell’s complete diagonal, \(\rho_{\mathrm{bnd}}\alpha_v(\Omega_P/{}^{(u)}\!a_P) A_{\mathrm{bnd}}/\ell_{\mathrm{bnd}}\); Denner’s boundary coefficient is given above. A boundary coefficient need not be obtained by copying the interior two-cell average.
5.9.2. Cell-Velocity and Face-Flux Corrections with SIMPLEC#
Subtracting the provisional momentum row from a corrected row, with coefficients held fixed, gives
SIMPLE neglects the neighbor-correction sum. SIMPLEC approximates each \(\vec u_N'\) by \(\vec u_P'\). The coefficient multiplying \(-\nabla p_P'\) in the cell-velocity correction is then \(\Omega_P/[{}^{(u)}\!a_P^{\mathrm{pred}}- \sum_N{}^{(u)}\!a_{P,N}]\). Inserting this cell-velocity correction into the selected face relation gives the corresponding face-flux correction.
For the one-step Pascau relation, with all history values held fixed but retaining the induced cell-velocity and cell-gradient corrections, subtraction gives
With the local-gradient SIMPLE cell-velocity correction, the first two terms cancel. SIMPLEC leaves an additional cell-gradient contribution. For the source-free interior continuity-subtracted row, \(\sum_N{}^{(u)}\!a_{P,N}=\widetilde a_P\) and \({}^{(u)}\!a_P^{\mathrm{pred}}=\widetilde a_P(1+\delta_P)/\alpha_v\). The ratio of the SIMPLEC and SIMPLE cell-velocity correction coefficients is therefore
Thus SIMPLEC multiplies the local SIMPLE cell-velocity correction by \(1+\kappa_P\), with
The remaining term in Eq. (5.9.2.2) is then \(-\alpha_v A_f\mathcal I_f[(\Omega/\widetilde a)\kappa\nabla p'] \cdot\hat{\vec n}_f\). To obtain a two-cell face formula, approximate this remaining cell-gradient term by the direct pressure jump, using \(\kappa_f=\alpha_v/(1-\alpha_v+\delta_f)\), giving the effective coefficient \((1+\kappa_f)d_f\). This is an additional face approximation. Source diagonals, boundaries, and temporal normalization require the corresponding coefficient balance. In the pressure-correction and algorithm chapters, \(d_f\) denotes the effective mass-flux correction coefficient for the selected algorithm, including this multiplier when SIMPLEC is used. The full algorithm comparison is developed in Chapter 7.
The selected momentum-interpolation family determines:
the provisional flux \(\dot m_f^{*,k+1}\),
the mass-flux correction coefficient \(d_f\), and
the explicit iteration, time-history, geometry, and source-balance corrections contained in \(\dot m_f^{*,k+1}\).
Once Eq. (5.9.2) has been obtained, each interior face adds the same diagonal/neighbor pattern to the pressure-correction matrix, regardless of the interpolation family. Chapter 6 develops that assembly, while Pressure-Correction Equation Forms describes the additional normalizations and compressible forms applied to the base coefficient.
5.10. Moving Grids and Relative Flux#
On a moving control volume, continuity uses the fluid flux relative to the grid:
where \(\dot{\Omega}_{g,f}\) is the swept grid-volume flux. Momentum interpolation constructs the fluid flux and its dependence on the adjacent-cell pressure difference. Continuity and its flux correction use the relative flux. At the predictor stage, the relative flux supplied to the pressure-correction equation is therefore
The geometric conservation law also requires the same time-consistent face area representation in:
the transported face flux,
the pressure and source terms in the momentum-interpolation formula, and
the swept grid-volume flux.
Otherwise a uniform flow may acquire a nonzero discrete residual solely from mesh motion. The space–time area construction used for this purpose is developed in Geometric Conservation Law for Moving Meshes.
5.11. Comparing the Completed Constructions#
For the fixed-grid, one-step scope of the derivations, the following table summarizes the pieces that must be used together. The pressure coefficients are the base SIMPLE mass-flux correction coefficients before further algorithm or continuity normalization.
Selection |
Provisional velocity contribution |
History or acceleration treatment |
Construction of \(d_f\) |
|---|---|---|---|
|
Ordinary geometric interpolation |
Previous-iteration face–cell flux difference |
Complete momentum diagonal |
|
Same as |
Classical iteration history plus an old-time face–cell difference with a common ratio |
Same direct coefficient as |
|
Interpolate \((1+\delta)\vec u^*\), then divide by \(1+\delta_f\) |
Weighted cell and face histories at both iteration and physical-time levels |
Non-temporal diagonal and \(1+\delta_f\) |
|
Same as |
Regrouped iteration term and separately reconstructed \(\vec T^k\) |
Same direct coefficient as |
|
Interpolate \([1+(1-\alpha_v)\delta]\vec u^*\), then divide by its face factor |
Weighted iteration history and separately reconstructed \(\vec T^*\) |
Non-temporal diagonal and \(1+(1-\alpha_v)\delta_f\) |
|
Projected-distance interpolation |
Iteration history and \(G_f\)-weighted old-time history |
\(\rho_f\alpha_v\beta_fF_f A_f/d_{PN}\), with matched density-weighted pressure/source residuals |
The opening selector table links each name to its derivation. Formulation names identify these coupled choices; comparing only their pressure coefficients omits part of the method.
5.12. Consistency Requirements#
A momentum-interpolation construction is consistent with the surrounding pressure-based method when:
the face flux responds directly to the pressure difference between the two cells sharing the face;
the face pressure-difference term and the cell pressure term removed from the interpolated velocities use the same coefficient;
nonlinear relaxation history is represented at the face as well as at the cells;
physical-time history is discretized consistently between cell velocity and face flux;
pressure, gravity, and other momentum sources use compatible discrete force and interpolation operators;
the base pressure-correction coefficient is the negative partial derivative of the selected mass flux with respect to the pressure difference, with the reconstructed cell quantities held fixed;
boundary coefficients respect whether pressure or flux is prescribed; and
moving-grid fluxes and area vectors satisfy the same geometric conservation construction.
The face flux carries the pressure, source, and history terms from the momentum equation into continuity. Its mass-flux correction coefficient supplies the local relation between a pressure-difference change and a flux change used to assemble the pressure-correction equation.
5.13. References#
Paul Bartholomew, Fabian Denner, Mohd Hazmil Abdol-Azis, Andrew Marquis, and Berend G. M. van Wachem. Unified formulation of the momentum-weighted interpolation for collocated variable arrangements. Journal of Computational Physics, 375:177–208, 2018. doi:10.1016/j.jcp.2018.08.030.
Seok-Ki Choi. Note on the use of momentum interpolation method for unsteady flows. Numerical Heat Transfer, Part A: Applications, 36(5):545–550, 1999. doi:10.1080/104077899274679.
J.H. Ferziger and M. Peric. Computational Methods for Fluid Dynamics. Springer, 3rd edition, 2002.
S. Majumdar. Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids. Numerical Heat Transfer, 13(1):125–132, 1988. doi:10.1080/10407788808913607.
A. Pascau. Cell face velocity alternatives in a structured colocated grid for the unsteady navier–stokes equations. International Journal for Numerical Methods in Fluids, 65:812–833, 2011. doi:10.1002/fld.2215.
C.L. Rhie and W.L. Chow. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation. AIAA Journal, 21:1525–1532, 1983.