Projection-Based Navier-Stokes Solvers

• Projection-Based Navier-Stokes solvers are algorithms that decouple velocity updates and pressure corrections via operator splitting to enforce divergence-free flow fields.
• They integrate advanced boundary conditions, high-order time-stepping, and stabilization techniques to efficiently tackle multiphase, multispace, and multiphysics fluid problems.
• Robust mathematical foundations, including error analysis and convergence studies, validate their use in scalable simulations of turbulence, fluid-structure interaction, and complex flow scenarios.

Projection-based Navier-Stokes solvers are a foundational class of algorithms for numerically integrating the incompressible Navier-Stokes equations, leveraging operator splitting between velocity update and pressure correction to enforce the divergence-free constraint at each time step. These schemes, originating with the Chorin–Temam approach, have evolved to incorporate advanced boundary condition treatments, stabilization protocols, high-order time-stepping, robust handling of multispace discretizations, and adaptation to multicomponent flows. Their capacity to decouple velocity and pressure enables efficiency, flexibility, and scalability, supporting applications ranging from canonical turbulence simulations to multiphysics fluid-structure and multiphase modeling.

1. Fundamental Time-Splitting Algorithms

The classical projection algorithm proceeds in three principal steps for each time increment $\Delta t$ (Pekker et al., 2024). Given a divergence-free velocity $u^n$, one computes:

• Momentum (Predictor) Step: Advance $u^n$ by solving

$$u^* = u^n + \Delta t \left[ - (u^n \cdot \nabla)u^n + \nu \Delta u^n + f^n \right]$$

applying appropriate boundary conditions to $u^*$.

• Pressure (Poisson) Correction: To enforce incompressibility, solve for $p^{n+1}$ such that

$$\nabla^2 p^{n+1} = \frac{1}{\Delta t}\nabla\cdot u^*$$

with necessary boundary conditions for $p^{n+1}$.

• Velocity Correction: Retrieve the divergence-free velocity

$u^{n+1} = u^* - \Delta t \nabla p^{n+1}$

resulting in $\nabla \cdot u^{n+1} = 0$.

The centered-in-time (Crank–Nicolson) pressure/convection variants can achieve second-order accuracy in time (Pekker et al., 2024). Extensions for higher-order implicit multistep schemes (BDF2, BDF3, etc.) are standard in level-2 codes (Dokken et al., 2019, Fehn et al., 2020).

2. Pressure–Tangential-Velocity Boundary Conditions

Recent analysis by Pekker & Pekker establishes a rigorously second-order-in-time projection framework for incompressible Navier-Stokes problems posed with prescribed pressure and tangential velocity at boundaries (open/free surfaces, inlets/outlets) (Pekker et al., 2024). The specification is:

• Dirichlet for boundary pressure: $p|_{\partial\Omega} = p_b(y, t)$
• Dirichlet for tangential velocity: $(u \cdot \tau)|_{\partial\Omega} = u_\tau^b(y, t)$
• Recovery of normal-velocity gradient: $$\partial_n (u \cdot n)|_{\partial\Omega} = -\partial_\tau (u \cdot \tau)|_{\partial\Omega}$$ (from $\nabla\cdot u = 0$)

With a centered pressure update (e.g., $\nabla P^{n+1/2}$), the boundary treatment preserves second-order temporal accuracy. If using lagged or forward pressure (fully explicit or implicit), accuracy is only first order at boundaries. This distinction is quantitatively demonstrated by normal-mode analysis of reflection coefficients. Implementation on staggered grids (MAC) simply modulates the Poisson matrix rows to enforce the Dirichlet pressure, with tangential velocity imposed in the predictor and incompressibility recovering the normal component (Pekker et al., 2024).

3. Multimesh, High-Order, and Discontinuous Galerkin Extensions

Projection-based solvers have been generalized to:

• Multimesh Finite Element Formulations: Multi-domain or overlap-union meshes are coupled via Nitsche-type weak enforcement of velocity continuity, with ghost-penalty stabilization ensuring coercivity and inf-sup stability. The time integration uses incremental pressure-correction splitting (BDF2-AB2). Convergence analysis demonstrates optimal rates, and performance is validated on Taylor-Green and Turek–Schäfer benchmarks (Dokken et al., 2019).
• High-Order Discontinuous Galerkin (DG) Discretizations: In DG schemes, the pressure gradient and velocity divergence coupling terms must be integrated by parts to avoid instabilities at coarse mesh or small $\Delta t$ regimes. Consistent boundary conditions for intermediate velocity and pressure correction (mirror principle) are essential. Mixed-order polynomial pairs (velocity degree $k$, pressure $k-1$) guarantee robust inf-sup stability; equal-order pairs require additional grad-div stabilization (Fehn et al., 2017). These practices restore projection-splitting to parity with monolithic DG in terms of accuracy and robustness.
• Arbitrary Lagrangian-Eulerian (ALE) Moving Meshes: Time-splitting and pressure-correction projection methods have been implemented in ALE-DG frameworks, preserving the geometric conservation law (GCL) by evaluating all weak forms and mass matrices on the mesh at the updated time level. Consistent boundary treatments for each projection substep are necessary for high-order accuracy (Fehn et al., 2020).

4. Projection Methods for Multiphase and Multiphysics Flows

Projection algorithms naturally extend to coupled fluid–interface problems:

• Cahn–Hilliard–Navier–Stokes (CHNS) Systems: Projection-based decoupling is used for phase-field multiphase flows, enforcing variable-density/viscosity corrections and consistent hydrodynamic–chemical-potential coupling in the velocity prediction and pressure-correction steps. The pressure-Poisson solver accommodates variable coefficients to handle density contrasts; mass conservation is ensured both globally and locally (Manna et al., 27 Aug 2025, Rabeh et al., 2024, Khanwale et al., 2021).
• Mass-Preserving Projection: For level-set two-phase simulations, projection steps are augmented with a global pressure mode found via a block Schur complement to enforce exact volume conservation, especially through morphological changes (breakup, coalescence) (Salac, 2015).
• Stochastic/Fluctuating Hydrodynamics: Projection splitting is adapted for thermal fluctuations via stochastic stress addition, requiring iterative correction steps for strict covariance matching at no-slip boundaries (Mancini et al., 2021).

5. Algorithmic and Discretization Strategies

A summary of distinctive discretization and implementation features:

Strategy	Key Attributes	Reference
Staggered/MAC	Pressure at centers, velocities at faces; robust divergence-free and boundary treatments.	(Pekker et al., 2024)
Collocated octree	Supra-convergent finite-difference stencils, sharp boundary interpolations, adaptive mesh refinement.	(Blomquist et al., 2023)
Multimesh FEM	Non-matching mesh coupling via Nitsche and ghost penalties, operator-split time stepping.	(Dokken et al., 2019)
DG and ALE	Integration-by-parts for coupling terms; mirror principle for boundary fluxes; high-order implicit/extrapolated time stepping; GCL preservation.	(Fehn et al., 2017, Fehn et al., 2020)
Bi-grid FEM	Implicit coarse-grid solve + stabilized high-mode semi-implicit correction on fine grid allows increased $\Delta t$ without loss of accuracy.	(Abboud et al., 2018)
Reduced-Order SBM	POD-Galerkin projection of snapshot-based solution spaces for rapid online evaluation; pressure–velocity inf-sup stabilization via stabilized full-order model.	(Karatzas et al., 2019)

6. Convergence, Stability, and Error Analysis

Mathematical foundations of projection methods are well-established:

• Existence and Stability: Discrete energy inequalities guarantee unconditional stability in normed spaces for both divergence-free projected velocity and intermediate velocity (Weber, 19 Aug 2025, Maeda et al., 2020, Abboud et al., 2018). Compactness lemmas confirm convergence of discrete sequences to Leray-Hopf weak solutions without elevated regularity assumptions (Weber, 19 Aug 2025, Maeda et al., 2020).
• Order of Accuracy: Temporal accuracy at boundaries hinges upon the alignment of BC specification and pressure update (centered splits yield $O(\Delta t^2)$, explicit/implicit splits reduce to $O(\Delta t)$) (Pekker et al., 2024). Spatial convergence rates align with the underlying finite element or difference schemes (e.g., $O(h^{k+1})$ for optimized FE pairs, second-order for node-based octree methods (Blomquist et al., 2023), fourth-order for high-order FV and method-of-lines approaches (Li et al., 2024)).
• Inf-sup Constraints and Pressure Modes: Discrete selection of velocity and pressure spaces (mixed order, penalty stabilization) is fundamental to the suppression of spurious pressure oscillations and to retaining optimal convergence (Fehn et al., 2017).

7. Practical Recommendations and Limitations

Best practices established by recent research include:

• Use centered-pressure ($\nabla P^{n+1/2}$) splits for maximizing time accuracy at boundaries with Dirichlet pressure–tangential velocity conditions.
• Carefully reconstruct intermediate velocities for consistent boundary enforcement and correct velocity divergence at the boundary; compute normal velocity corrections via the incompressibility constraint (Pekker et al., 2024).
• For discontinuous Galerkin and ALE variants, employ integration-by-parts of divergence and gradient operators, mirror principle for BCs, and mixed-order element pairs (Fehn et al., 2017, Fehn et al., 2020).
• In bi-grid, stabilized, or monolithic meshfree approaches, apply high-mode damping or overdetermined least-squares to ensure local mass conservation and algorithmic robustness (Abboud et al., 2018, Suchde et al., 2017).
• For multiphase and variable-density flows, ensure explicit coupling in pressure-Poisson equations and enforce global mass conservation via augmented pressure modes if needed (Salac, 2015).
• In reduced-order modeling, propagate inf-sup stability from the full-order stabilized solver to the reduced subspace; pressure supremizer enrichment may be omitted if the base scheme internalizes stabilization (Karatzas et al., 2019).

These projection-based methods are implemented in major open-source CFD packages (e.g., Gerris/Basilisk), and have shown performance advantages over fully coupled saddle-point solvers in large-scale, adaptive, and multiphase applications.

---

References: (Pekker et al., 2024, Dokken et al., 2019, Fehn et al., 2017, Blomquist et al., 2023, Abboud et al., 2018, Weber, 19 Aug 2025, Maeda et al., 2020, Fehn et al., 2020, Salac, 2015, Khanwale et al., 2021, Manna et al., 27 Aug 2025, Rabeh et al., 2024, Shirokoff et al., 2010, Karatzas et al., 2019, Azaïez et al., 11 Dec 2025, Li et al., 2024, Mancini et al., 2021, Suchde et al., 2017).