Proposed Changes

The current work (part of GSoC) provides a working version of a pressure-based algorithm for the incompressible flow solver as an alternative to the existing Density-based solver. Below, the reader may find the algorithm which has been implemented, as well as the current progress of the code and challenges. All the way at the bottom one can find performance comparisons between the DB and PB solvers for some test cases.

Algorithm

A lot of versions of pressure-based algorithms exist, and many different versions can be implemented. Here, we opt for versions of the original SIMPLE/PISO algorithm, it is briefly defined here for clarity.

First, the momentum equations are solved, starting from the previous time step's momentum $\rho \vec{u}^{(0)}$, pressure $p^{(0)}$, and face flux $\rho \vec{u}_f^{(0)}$. The resulting momentum is the predicted momentum, here its discretized form is shown, as its coefficients are used in the subsequent equations

$A_p(\rho \vec{u})_p^{(1)}+\sum_n A_n (\rho \vec{u})^{(1)}_n=-V\nabla p^{(0)}+S_m$

The subsequent momentum is not necessarily incompressible, the pressure correction equation can be derived by rewriting it as follows, using a term often called H by A

$(\rho \vec{u})_p^{(1)}=-\frac{\sum_n A_n (\rho \vec{u})_n^{(1)}}{A_p}-\frac{V}{A_p}\nabla p+S_m=\frac{H\left(\rho \vec{u}^{(1)}\right)}{A}-\frac{V}{A_p}\nabla p^{(0)}+S_m$

$(\rho \vec{u})_p^{(2)}=\frac{H\left(\rho \vec{u}^{(1)}\right)}{A}-\frac{V}{A_p}\nabla p^{(1)}+S_m$

Note how a simplification is used here where the HbyA term is neglected. Subtracting these two equations yields the first pressure correction equation for $p'$ as

$\nabla \cdot \left( \frac{V}{A_p}\nabla p'\right)=\nabla \cdot (\rho \vec{u})^{(1)}$

Make note that for the divergence here, we require the face mass fluxes, which are computed using Rhie-Chow interpolation to avoid odd-even decoupling. After the equation is solved, using the pressure correction $p'$, the pressure and momentum are corrected according to

$p^{(1)} = p^{(0)} + p' ,\quad \rho u^{(2)} = \rho u^{(1)}+\rho u'.\quad \rho u'= -\frac{V}{A_p}\nabla p'$

So far, this is equal to a pseudo-transient version of the SIMPLE algorithm. This algorithm however suffers from a very tight stability condition on the time-step size. Therefore, multiple pressure corrections can be applied, which for two corrections is originally called the PISO algorithm.

The second pressure correction does not neglect the HbyA term, which then results in the equation

$\nabla \cdot \left( \frac{V}{A_p}\nabla p'\right)=\nabla \cdot (\rho \vec{u})^{(2)}+\nabla\cdot\left(\frac{H\left(\rho \vec{u}'\right)}{A}\right)$

And the new correction equations are defined as

$p^{(2)} = p^{(1)} + p' ,\quad \rho u^{(3)} = \rho u^{(2)}+\rho u'.\quad \rho u'= \frac{H\left(\rho \vec{u}'\right)}{A}-\frac{V}{A_p}\nabla p'$

Note that HbyA here uses the previous momentum correction and is thus the same quantity as the one used in the second pressure correction equation. Later pressure correction equations follow analogously.

Progress:

• Pressure-based solver added as alternative to density-based solver for the incompressible flow equations.
• The pressure-based solver has only been tested for constant density cases for basic Navier-Stokes and Euler flow.
• The pressure-based solver is implemented based on a pseudo time-stepping approach to remain consistent with the other solvers in SU2. The pressure-based solver is currently set to the SIMPLE algorithm by default, with options for SIMPLEC and PISO available.
• The Poisson solver is a major bottleneck in the computation speed. To account for this, an option was added to use a different linear solver and preconditioner for the poisson solver.
• For details on the implementation of the algorithm and the responsibility distribution please see the file CPBFluidIteration.cpp.

Issues

Performance:

• Convergence issues with RANS (SA and SST) on fine meshes with high Reynolds numbers.
• The Poisson solver can sometimes struggle a lot due to high Reynolds numbers and fine meshes, and thus require a ridiculous number of iterations to converge reasonably. Possible fixes include adding multigrid support or a DIC preconditioner (far less efficient). Multigrid support is tricky as SU2 currently only considers multigrid for the main (flow) solver and not for auxiliary solvers.
• Periodic boundary conditions have not been implemented/tested at all as of yet.
• Any code related to adjoints has not been considered at all either.

Code:

• Parallelization with OMP gives wrong results, MPI however does work as expected.

TODO list

• Fix issues mentioned above
• Docs page

Related Work

This work is based on earlier attempts by Nitish Anand (2024) and Akshay Koodly (2021), see feature branches feature_PBFlow_V8 and feature_Pressure_based respectively. Also see PR #2210

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

Result showcase

Convergence history of the inviscid flow around a hydrofoil at a 5 degree aoa.

The pressure coefficient along the surface of the hydrofoil at a 5 degree aoa and the corresponding lift coefficients, X-FOIL predicts C_L=0.6.

Convergence history of the lid driven cavity problem, note that CFL=60 is the highest stable CFL for the PB solver, whereas the DB solver does not have this CFL related stability issue.

The skin friction coefficient for turbulent flow over a (rough) flat plate with SA.