Motivation

For steady-state incompressible flow problems, classical time-marching schemes use a single global time step \(\Delta t\) shared by all computational cells. This global step is limited by the most restrictive cell (typically a small cell or a high-velocity region), leading to very slow convergence toward the steady state.

The Implicit_Euler_Steady_Scheme addresses this by replacing the global time step with a local (per-face) pseudo-time step \(\Delta\tau_f\), based on the local convective CFL condition at each face \(f\):

\[\Delta\tau_f = c_s \cdot f \cdot \frac{V_f}{|u_f|} \]

where \(V_f\) is the face-associated volume, \(|u_f|\) the convective velocity magnitude, \(c_s\in(0,1]\) a CFL safety factor (default \(c_s=0.8\)) and \(f\geq1\) a dynamic amplification factor (keyword facsec, default 1, bounded by facsec_max). Each cell can therefore advance at its own optimal pseudo-time step, dramatically accelerating convergence on heterogeneous meshes.

Mathematical formulation

The incompressible Navier-Stokes equations in discrete form read:

\[\begin{aligned} \mathbf{M}\frac{\mathbf{U}^{n+1}-\mathbf{U}^n}{\Delta\tau}+\mathbf{A}\mathbf{U}^{n+1}-\mathbf{L}(\mathbf{U}^n)\mathbf{U}^{n+1}+\mathbf{B}^T P^{n+1} &= \mathbf{S}\\ \mathbf{B}\,\mathbf{U}^{n+1} &= 0 \end{aligned} \]

where \(\mathbf{M}\) is the mass matrix, \(\mathbf{A}\) the diffusion operator, \(\mathbf{L}(\mathbf{U}^n)\) the linearised convection operator evaluated at \(\mathbf{U}^n\), \(\mathbf{B}\) and \(\mathbf{B}^T\) the discrete divergence and gradient operators, \(P\) the pressure and \(\mathbf{S}\) the source term. In the local-time-step formulation, the scalar \(\Delta\tau\) is replaced by the diagonal matrix \(\boldsymbol{\Delta\tau}=\text{diag}(\Delta\tau_f)\), giving the modified inertia term \(\mathbf{M}\,\boldsymbol{\Delta\tau}^{-1}\).

The local time step at face \(f\) is \(\Delta\tau_f=c_s\min_{\text{eq}}(\Delta\tau_f^{\text{CFL, eq}})\), the minimum being taken over all equations of the problem (e.g. Navier-Stokes, \(k\)- \(\varepsilon\)) to ensure simultaneous stability; only the convective CFL constraint is used (diffusion is treated implicitly). The stationarity criterion is the relative change in the velocity field between two consecutive pseudo-time steps, \(\frac{\|\mathbf{U}^{n+1}-\mathbf{U}^n\|}{\|\mathbf{U}^n\|}<\varepsilon_{\text{statio}}\), which is independent of the magnitude of the local pseudo-time steps.

Predictor-corrector algorithm

The solver Implicit_steady implements a projection method (SIMPLE-like) with local time steps, derived from the Simpler solver.

Predictor (velocity). Assemble the momentum system with inertia \(\mathbf{M}\,\boldsymbol{\Delta\tau}^{-1}\) and solve for the intermediate velocity \(\mathbf{U}^*\):
\[\left(\frac{\mathbf{M}}{\boldsymbol{\Delta\tau}}+\mathbf{A}-\mathbf{L}(\mathbf{U}^n)\right)\mathbf{U}^*=\frac{\mathbf{M}}{\boldsymbol{\Delta\tau}}\mathbf{U}^n-\mathbf{B}^T P^n+\mathbf{S} \]
Pressure correction. Compute the divergence of \(\mathbf{U}^*\) and assemble the pressure equation:
\[\left(\mathbf{B}\,(\mathbf{M}/\boldsymbol{\Delta\tau})^{-1}\,\mathbf{B}^T\right)P'=\mathbf{B}\,\mathbf{U}^* \]
Velocity correction and pressure update.
\[\mathbf{U}^{n+1}=\mathbf{U}^*-(\mathbf{M}/\boldsymbol{\Delta\tau})^{-1}\mathbf{B}^T P', \qquad P^{n+1}=P^n+P' \]

The pressure matrix is reassembled only when both conditions hold: (i) the maximum relative variation of \(\boldsymbol{\Delta\tau}\) exceeds 25%, and (ii) the time-step index is a multiple of the reassembly interval (default every 100 steps). The 25% threshold is the safe upper bound derived from \(c_s=0.8\) (a variation up to \(1/c_s-1=25\%\) can be absorbed without violating the CFL constraint). When neither condition is met, \(\boldsymbol{\Delta\tau}\) and the pressure matrix are kept from the previous iteration, ensuring all three steps use the same \(\boldsymbol{\Delta\tau}\).

User parameters

Keyword	Type	Default	Description
seuil_statio	double	\(10^{-12}\)	Stationarity threshold
nb_pas_dt_max	int	\(\infty\)	Maximum number of pseudo-time steps (recommended for RANS cases — see note)
facsec_max	double	1	Upper bound for the dynamic amplification factor \(f\). With the default 1, \(f\) stays fixed and the local time step is purely CFL-based; when > 1, \(f\) grows automatically (based on residual reduction) up to facsec_max.

Note: RANS turbulence models. With \(k\)- \(\varepsilon\) or \(k\)- \(\omega\) models, the turbulence variables may exhibit persistent low-amplitude oscillations that prevent the stationarity criterion from being satisfied, causing the simulation to run indefinitely. In such cases, set nb_pas_dt_max as the primary stopping criterion and monitor convergence through the printed residuals.

Data file syntax

Implicit_Euler_Steady_Scheme sch
Read sch {
    nb_pas_dt_max  500
    seuil_statio   1.e-6
    facsec_max     10    # optional: enables dynamic facsec growth up to 10
    Solveur Implicit_steady {
        solveur_pression GCP { precond ssor { omega 1.5 } seuil 1.e-10 }
    }
}