Segregated Flow Solver

Discretization of the Momentum Equation

The convective and diffusive fluxes are discretized as described in Convective Flux and Diffusive Flux. The discretization of the pressure gradient requires special attention.

Face Pressure Evaluation

Pressure appears as pressure gradient in the momentum equations. In order to compute the pressure gradient term in Eqn. (922), the pressure must be evaluated at each face.

where $\overline{a_0}$ and $\overline{a_1}$ are the average of the momentum coefficients for all components of momentum for cells 0 and 1, respectively. $p_{f0}$ and $p_{f1}$ are interpolated from cell values and reconstruction gradients according to Eqn. (911).

When the delta-V dissipation expert property is turned on, a dissipation term is incorporated into the face pressure evaluation, as follows:

where ${\mathbf{\text{v}}}_{f0}$ and ${\mathbf{\text{v}}}_{f1}$ are also evaluated according to Eqn. (911).

At boundaries, the cell 1 contribution comes from a fictitious “ghost cell” in which the velocity has been reflected about the boundary face. Therefore Eqn. (923) becomes:

With delta-V dissipation on, Eqn. (924) becomes:

at a boundary.

Pressure-Correction Equation

The pressure-velocity coupling is achieved by re-writing the continuity equation in terms of a mass flux correction ${{\dot{m}}^{\prime }}_f$:

The uncorrected face mass flux ${\dot{m}}_f^{\text{*}}$ is computed after the discrete momentum equations have been solved. The momentum equations are solved with a guessed pressure field $p^{*}$, which does not satify continuity. The mass flow correction ${{\dot{m}}^{\prime }}_f$ is required to satisfy continuity.

The uncorrected mass flux at an interior face can be written in terms of the cell variables as follows:

where ${\mathbf{\text{v}}}_0^{\text{*}}$ and ${\mathbf{\text{v}}}_1^{\text{*}}$ are the cell velocities after the discrete momentum equations have been solved. Linear interpolation of the velocities at the cell centers to the face results in unphysical checkerboarding of pressure. To prevent this checkerboarding, the Rhie-and-Chow dissipation${\Upsilon }_f$ is introduced:

where:

In multiphase flows, due to presence of interaction between the phases, the momentum coefficients form a matrix. However, only the diagonal terms of the matrix are considered. These are represented by ${\overline{a}}_0$ and ${\overline{a}}_1$ in Eqn. (930).

Unsteady Dissipation Correction

For either coefficient $\overline{a}={\overline{a}}_0\text{or}{\overline{a}}_1$ , when unsteady dissipation corrections are enabled, the coefficients $\overline{a}$ are modified based on an Acoustic CFL.

• Let Acoustic CFL be ${CFL}_a=c\frac{\Delta t}{\Delta x}$ where $c$ is the speed of sound, $\Delta t$ is the simulation time-step, and $\Delta x$ is the local grid size.
• Let $\text{limCFL}$ be the Limit Acoustic CFL value specified by the user.

Then:

• If $\text{CFL}>\text{limCFL}$, $\overline{a}$ is used without modification.
• Otherwise, $\overline{a}$ is replaced with $\overline{a}{}^{\prime }$, as follows:

  $\overline{a}{}^{\prime }={\overline{a}}_u(1+\frac{{\text{CFL}}_r{\overline{a}}_{uns}}{\text{max}({\overline{a}}_u,{1\times 10}^{-20})})$

  where:

  • ${\overline{a}}_u={\overline{a}-\overline{a}}_{uns}$ and
${\overline{a}}_{uns}=\frac{3}{2}\rho \frac{\text{Vol}}{\mathrm{\Delta }t}$
  • ${\text{CFL}}_r=\text{min}(1.0,\frac{\text{CFL}}{\text{limCFL}})$

$p_L$ and $p_R$ are the pressures values computed by linear reconstruction at two symmetric points $({\mathbf{\text{x}}}_L,{\mathbf{\text{x}}}_R)$ with regard to the face-centroid, in the direction of face normal.

${\mathbf{\text{x}}}_L={\mathbf{\text{x}}}_f-d_s\mathbf{\text{n}}$ and ${\mathbf{\text{x}}}_R={\mathbf{\text{x}}}_f+d_s\mathbf{\text{n}}$

where:

• $\mathbf{\text{n}}$ is the unit vector associated with area vector ${\mathbf{\text{a}}}_f$.
• ${\mathbf{\text{x}}}_f$ is the position vector of the face centroid.
• $d_s$ is a scalar distance, where:
  • If $\text{max}(\frac{V_1}{V_0},\frac{V_0}{V_1})>10$, then $d_s=\frac{0.5}{\left|\mathbf{\text{α}}\right|}$.

  • Otherwise, $d_s=\text{min}(d_{s0},d_{s1})$, where:

    • $d_{s0}=d{\mathbf{\text{x}}}_0\cdot \mathbf{\text{n}}$ (the projection of $d{\mathbf{\text{x}}}_0={\mathbf{\text{x}}}_f-{\mathbf{\text{x}}}_0$ onto the face area normal)
    • $d_{s1}=d{\mathbf{\text{x}}}_1\cdot \mathbf{\text{n}}$ (the projection of $d{\mathbf{\text{x}}}_1={\mathbf{\text{x}}}_f-{\mathbf{\text{x}}}_1$ onto the face area normal).

Density Correction

For compressible flow, the density correction $\mathrm{\rho \prime }$ is introduced:

$$\begin{array}{c}{\dot{m}}_f={(\rho +\mathrm{\rho \prime })}_f({\text{v}}_{fn}^{\text{*}}+{\text{v}\text{'}}_{fn})\left|\mathbf{\text{a}}\right|\\ =({\rho }_f{v}_{fn}^{\text{*}}+{\mathrm{\rho \prime }}_f{v}_{fn}^{\text{*}}+{\rho }_f{\mathrm{v\prime }}_{fn}+{\mathrm{\rho \text{'}}}_f{\mathrm{v\prime }}_{fn})\left|\mathbf{\text{a}}\right|\end{array}$$

(932)

where the subscript “fn” denotes the face-normal component.

$${\rho }_f{\text{v}\text{'}}_{fn}\left|\mathbf{\text{a}}\right|\equiv -Q_f\left({p\prime }_1-{p\prime }_0\right)$$

(933)

where ${p\prime }_1$ and ${p\prime }_0$ are the cell pressure corrections, and:

$${\mathrm{\rho \prime }}_f{\text{v}}_{fn}^{\text{*}}\left|\mathbf{\text{a}}\right|=\frac{{\dot{m}}_f^{\text{*}}}{{\rho }_f}{\left(\frac{\partial \rho }{\partial p}\right)}_T{\mathrm{p\prime }}_{\text{upwind}}$$

(934)

where:

$${\mathrm{p\prime }}_{\text{upwind}}=\{\begin{array}{c}{p\prime }_0\text{ for }{\dot{m}}_f^{\text{*}}>0\\ {p\prime }_1\text{ for }{\dot{m}}_f^{\text{*}}<0\end{array}$$

(935)

Since you can choose to neglect the second-order term ${\mathrm{\rho \text{'}}}_f{v}_{fn}\text{'}$, these equations can be combined to define the mass flow correction as:

$${{\dot{m}}^{\prime }}_f=Q_f\left({p^{\prime }}_0-{p^{\prime }}_1\right)+\frac{{\dot{m}}_f^{\text{*}}}{{\rho }_f}{\left(\frac{\partial \rho }{\partial p}\right)}_T{\mathrm{p\prime }}_{\text{upwind}}$$

(936)

Discrete Pressure Correction Equation

The discrete pressure correction equation is obtained from Eqn. (927) and Eqn. (936) and is written in coefficient form as:

$$p{\prime }_p+{\displaystyle \sum _n^{}}{a}_{n}p{\prime }_n=r$$

(937)

The residual r is the net mass flow into the cell:

$$r=-{\displaystyle \sum _n^{}}{\dot{m}}_f^{\text{*}}$$

(938)

SIMPLE Algorithm

The SIMPLE algorithm is used to control the overall solution. This algorithm is summarized as follows:

1. Set the boundary conditions.
2. Compute the reconstruction gradients of velocity and pressure.
3. Compute the velocity and pressure gradients.
4. Solve the discretized momentum equation. This creates the intermediate velocity field ${\mathbf{\text{v}}}^{\text{*}}$ .
5. Compute the uncorrected mass fluxes at faces ${\dot{m}}_f^{\text{*}}$ .
6. Solve the pressure correction equation. This produces cell values for the pressure correction $\mathrm{p\prime }$ .
7. Update the pressure field:

   $$p^{n+1}=p^n+\omega p\text{'}$$

   (939)

   where $\omega$ is the under-relaxation factor for pressure.

8. Update the boundary pressure corrections ${p\prime }_b$ .
9. Correct the face mass fluxes:

   $${\dot{m}}_f^{n+1}={\dot{m}}_f^{\text{*}}+{\dot{m}\prime }_f$$

   (940)
10. Correct the cell velocities:

    $$v_p^{n+1}=v_p^{*}-\frac{V\nabla \mathrm{p\prime }}{{\mathbf{\text{a}}}_p^v}$$

    (941)

    where $\nabla p\text{'}$ is the gradient of the pressure corrections, ${\mathbf{a}}_p^{\text{v}}$ is the vector of central coefficients for the discretized velocity equation, and $V$ is the cell volume.

11. Update density due to pressure changes.
12. Free all temporary storage.

SIMPLE Consistent (SIMPLEC) Algorithm

SIMPLEC follows the same steps as summarized for the SIMPLE algorithm. It differs in the form of the equation for the pressure correction $p^{\prime }$ . In SIMPLE, the velocity correction equation at a cell is written in terms of pressure correction, neglecting velocity corrections at the neighbouring cells. Taking divergence of this equation and making use of continuity gives the $p^{\prime }$ equation.

In its velocity correction equation, SIMPLEC assumes that the velocity corrections in the neighboring cells are equal to those at the central cell. This results in a different equation for $p^{\prime }$ .

Another difference is that the pressure under-relaxation relating to time-step size and velocity under-relaxation is implicitly built into the $p^{\prime }$ equation, so that there is no explicit pressure under-relaxation as in Eqn. (939).

In SIMPLEC, Eqn. (941) becomes:

where ${\mathbf{a}}_n^v$ are the velocity neighbor coefficients.

PISO Algorithm

The PISO algorithm is suited for transient cases where the convective Courant number is small.

Comparing PISO with SIMPLE:

• PISO is faster than SIMPLE at short time-steps, though both algorithms have the same level of temporal accuracy. [254]
• PISO becomes unstable at long time-steps, when the combined CFL rises much above 10, while SIMPLE remains stable.
• As time-step size increases, SIMPLE loses temporal accuracy of transient solutions. However SIMPLE can still obtain accurate steady state solutions, if they exist, by using large time-step size.

Therefore, use SIMPLE for cases where large time-step sizes can be used, such as in achieving steady state solution, but use PISO where the problem demands use of small time-step size and greater temporal accuracy.

The PISO algorithm shares many features with the SIMPLE algorithm. The algorithm is summarized as follows:

1. Set the boundary conditions.
2. Compute the reconstruction gradients of velocity and pressure.
3. Compute the velocity and pressure gradients.
4. Solve the discretized momentum equation. This action creates the intermediate velocity field $v^{*}$ .
5. Compute the uncorrected mass fluxes at faces ${\dot{m}}_f^{*}$ .
6. Construct the pressure correction equation and test for convergence of the PISO algorithm. If PISO has converged, then set the last iteration flag to true.
7. Solve the pressure correction equation. This action produces cell values for the pressure correction $p^{\prime }$ .
8. Update the pressure field:

   $$p^{n+1}=p^n+{\omega p}^{\prime }$$

   (943)

   where $\omega$ is the under-relaxation factor for pressure. If this iteration is the last, $\omega$ is equal to one.

9. Update the boundary pressure corrections ${p^{\prime }}_b$ .
10. Correct the face mass fluxes:

    $${\dot{m}}_f^{n+1}={\dot{m}}_f^{*}+{{\dot{m}}^{\prime }}_f$$

    (944)
11. Correct the cell velocities:

    $$v_p^{n+1}=v_p^{*}-\frac{V\omega \nabla p^{\prime }}{a_p^v}$$

    (945)

    where $\nabla p^{\prime }$ is the gradient of the pressure corrections, $a_p^v$ is the vector of central coefficients for the discretized linear system representing the velocity equation, and $V$ is the cell volume.

12. Update density due to pressure changes.
13. Solve any PISO-coupled equations such as for Species or Energy.
14. Recalculate density from the equation of state (except for the last corrector).
15. If this iteration is the last, then go to step 18.
16. Explicitly update the momentum equation with corrections due to the neighbor cells.
17. Go to step 5 and continue until convergence is reached.
18. Free all temporary storage.