Add Lid Driven Cavity Flow tutorial for pressure-based solver - #218
Open
thijsaalbers wants to merge 6 commits into
Open
thijsaalbers wants to merge 6 commits into
thijsaalbers wants to merge 6 commits into
Conversation
This was referenced Aug 23, 2026
pcarruscag
added a commit
to su2code/SU2
that referenced
this pull request
Sep 15, 2026
## 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 velocity $\vec{u}^{(0)}$, pressure $p^{(0)}$, and face
velocity $\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})_n^{(1)}=-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((\rho\vec{u})^{(1)})}{A}-\frac{V}{A_p}\nabla p^{(0)}+S_m$
$(\rho\vec{u})_p^{(2)}=\frac{H((\rho\vec{u})^{(1)})}{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 \left(
\rho \vec{u}^{(1)}\right)$
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 \left(
\rho \vec{u}^{(2)}\right)+\nabla\cdot\left(\frac{H((\rho
u)')}{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((\rho u)')}{A}-\frac{V}{A_p}\nabla p'$
Note that HbyA here uses the previous velocity 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.
- A tutorial showcasing the different options of the PB solver for the
lid driven cavity flow problem has been added, please see
su2code/Tutorials#86 and su2code/su2code.github.io#218
- Regression tests have been added.
### Future work:
#### Performance:
- 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.
- Convergence issues with RANS (SA and SST) on fine meshes with high
Reynolds numbers. Tests have shown that cases such as flow over a flat
plate converges fine. However, external aerodynamic cases such as the
naca0012 RANS test case do not converge well at all. The convergence
does slightly improve when we switch out the mesh for a more uniform
unstructured mesh without large aspect ratio cells in the wake of the
airfoil, although this only slightly helps. The flat plate turbulence
test case also uses large aspect ratio cells so this is not the sole
issue. The RANS solver also often requires many iterations of the
Poisson solver to converge reasonably, this is however not the reason
for the lack of convergence.
- 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 due to a unknown issue in
the Poisson solver, MPI however does work as expected.
- The restart solution currently does not write the edge mass fluxes
which are used by the pressure-based solver. Therefore restarted
solutions have to estimate these mass fluxes based on the average of the
nodal solutions. This results in the restarted solutions starting from
slightly different residuals. A fix for this issue can be to let the
edge mass fluxes be stored in the restart file.
- Support for the energy equation is added trough both the enthalpy
equation as well as trough the weakly coupled heat solver. This support
has however not yet been verified trough actual test cases.
- The numerics class for the convective residuals write the continuity
parts of the flux and jacobians. This will allow for a coupled solver in
the future, but as it stands now this code is _not_ used.
### TODO list
- Fix issues mentioned above (left for future work)
## 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
- [x] I am submitting my contribution to the develop branch.
- [x] My contribution generates no new compiler warnings (try with
--warnlevel=3 when using meson).
- [x] My contribution is commented and consistent with SU2 style
(https://su2code.github.io/docs_v7/Style-Guide/).
- [x] I used the pre-commit hook to prevent dirty commits and used
`pre-commit run --all` to format old commits.
- [x] I have added a test case that demonstrates my contribution, if
necessary.
- [x] I have updated appropriate documentation (Tutorials, Docs Page,
config_template.cpp), if necessary.
## Result showcase
### Inviscid Hydrofoil
<img
src="https://github.com/user-attachments/assets/2178b449-557d-4ec5-89b6-0d49ace6b612"
width="600">
*Convergence history of the inviscid flow around a hydrofoil at a 5
degree aoa.*
<img
src="https://github.com/user-attachments/assets/cd28de03-1c0a-4927-a5a3-567c454eb5c2"
width="600">
*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.*
### Lid Driven Cavity
<img
src="https://github.com/user-attachments/assets/0433d13c-b76d-4b9d-ac2c-615539877c9e"
width="600">
*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.*
### Flatplate RANS
<img
src="https://github.com/user-attachments/assets/b8db128d-e919-46a5-9b9d-aa17c1b94219"
width="600">
*The skin friction coefficient for turbulent flow over a (rough) flat
plate with SA.*
---------
Co-authored-by: Pedro Gomes <pcarruscag@gmail.com>
Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds a tutorial to showcase how the options of the pressure-based solver can be used. Please see the corresponding solver su2code/SU2#2812 which it is a part of.