Solver Configuration and Numerical Methods

Scope

This page documents the governing equations, discretisation schemes, solver algorithms, and physics models available in AortaCFD. Three pre-configured numerics profiles --- robust, standard, and precise --- allow practitioners to trade numerical diffusion against stability depending on mesh quality and simulation objectives. All settings are compatible with OpenFOAM 12 and the foamRun solver modules.

---

1. Governing Equations

AortaCFD solves the incompressible Navier--Stokes equations for blood flow in patient-specific aortic geometries.

1.1 Continuity (Mass Conservation)

\[ \nabla \cdot \mathbf{u} = 0 \]

1.2 Momentum Conservation

\[ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{u}) = -\frac{1}{\rho}\,\nabla p + \nabla \cdot \bigl(\nu_{\text{eff}}\,\nabla \mathbf{u}\bigr) \]

where:

Symbol	Description	Value / Note
\(\mathbf{u}\)	Velocity vector field	Solved variable
\(p\)	Pressure (kinematic or static)	Solved variable
\(\rho\)	Blood density	1060 kg m\(^{-3}\)
\(\nu_{\text{eff}}\)	Effective kinematic viscosity	\(\nu + \nu_t\) (RANS) or \(\nu + \nu_{\text{sgs}}\) (LES)
\(\nu\)	Molecular kinematic viscosity	\(3.7736 \times 10^{-6}\) m\(^{2}\) s\(^{-1}\)

Blood is modelled as an incompressible Newtonian fluid with constant properties (\(\rho = 1060\) kg m\(^{-3}\), \(\mu = 4.0 \times 10^{-3}\) Pa s), consistent with the approach widely adopted in aortic computational fluid dynamics Steinman 2002; Taylor & Figueroa 2009. For laminar flow, \(\nu_{\text{eff}} = \nu\); for Reynolds-Averaged Navier--Stokes (RANS) simulations, \(\nu_{\text{eff}} = \nu + \nu_t\), where \(\nu_t\) is the turbulent eddy viscosity.

---

2. Spatial Discretisation

AortaCFD employs the finite volume method (FVM) Ferziger & Perić 2002 with a collocated variable arrangement on unstructured polyhedral meshes generated by snappyHexMesh. The integral form of the governing equations is discretised over each control volume, yielding a system of algebraic equations.

2.1 Gradient Schemes

Gradient reconstruction is performed via the Gauss theorem with optional cell-based limiting to prevent unbounded overshoots on highly skewed cells.

Scheme	Formal Order	Behaviour	Used In
Gauss linear	2nd	Unbounded, exact on uniform meshes	---
cellLimited Gauss linear 1.0	2nd	Full Barth--Jespersen limiting	robust
cellLimited Gauss linear 1	2nd	Moderate limiting	standard, precise

The limiter coefficient ranges from 0 (no limiting, pure second order) to 1 (full limiting, maximum boundedness). The robust profile applies full limiting to prevent gradient overshoots on cells with high skewness.

2.2 Convection (Divergence) Schemes

The convective term \(\nabla \cdot (\mathbf{u} \otimes \mathbf{u})\) is the primary source of numerical diffusion and potential instability. AortaCFD provides three convection schemes, selected automatically by profile:

Scheme	Order	Properties	Profile
Gauss upwind	1st	Bounded, monotone, diffusive	robust
Gauss linearUpwind default	2nd	Low diffusion, bounded	standard
Gauss LUST grad(U)	2nd	Linear-Upwind Stabilised Transport	precise

Total Variation Diminishing (TVD) Limiters

The limitedLinearV and limitedLinear schemes apply TVD limiting Jasak 1996 to scalar and vector fields, respectively. TVD schemes guarantee that no new local extrema are introduced during the convective update, preserving boundedness while maintaining second-order accuracy in smooth regions Jasak 1996; Versteeg & Malalasekera 2007. The coefficient of 1 enforces full limiting for maximum monotonicity.

The LUST (Linear-Upwind Stabilised Transport) scheme blends 75% central differencing (Gauss linear) with 25% linear-upwind stabilisation. This yields substantially lower numerical diffusion than pure linearUpwind, making it suitable for Large Eddy Simulation where resolved turbulent structures must not be artificially damped.

2.3 Laplacian Schemes

The diffusive (Laplacian) term discretisation must account for mesh non-orthogonality. On a perfectly orthogonal mesh the standard two-point stencil is exact; on non-orthogonal meshes an explicit correction is required Ferziger & Perić 2002; OpenFOAM 2024.

Scheme	Correction	Profile
Gauss linear limited corrected 0.5	Limited non-orthogonal correction	robust, standard
Gauss linear limited corrected 0.5	Tuned limited correction	precise

The limited corrected formulation blends between the uncorrected and fully corrected Laplacian. A coefficient of 0.5 provides stability on meshes with non-orthogonality up to approximately 70 degrees; higher coefficients (closer to 1.0) provide more accurate correction but may introduce instability on poor-quality meshes OpenFOAM 2024.

2.4 Surface-Normal Gradient Schemes

Surface-normal gradients (snGradSchemes) are set consistently with the Laplacian correction:

Profile	Scheme
robust	limited corrected 0.5
standard	limited 0.5
precise	limited 0.5

---

3. Temporal Discretisation

Scheme	Order	Stability	Profile
Euler	1st	Unconditionally stable (A-stable)	robust
backward	2nd	Unconditionally stable, fully implicit	standard
CrankNicolson 0.9	2nd	Nearly centred; off-centring parameter 0.9 for stability	precise

Crank--Nicolson Off-Centring

Pure Crank--Nicolson (\(\alpha = 1.0\)) is second-order accurate but can exhibit oscillatory behaviour on coarse temporal grids. The coefficient \(\alpha = 0.9\) blends 90% Crank--Nicolson with 10% Euler implicit, retaining second-order formal accuracy while improving robustness. This provides superior phase accuracy compared to the fully implicit backward scheme, which is important for resolving pulse-wave propagation in aortic haemodynamics.

All profiles employ adaptive time-stepping controlled by the maximum Courant number:

Profile	Max Co	Initial \(\Delta t\)	Max \(\Delta t\)
robust	1.0	\(10^{-6}\) s	\(10^{-3}\) s
standard	0.8	\(10^{-6}\) s	\(10^{-3}\) s
precise	0.5	\(10^{-6}\) s	\(8 \times 10^{-4}\) s

The conservative initial time step (\(10^{-6}\) s) prevents Courant-number spikes during flow startup.

---

4. PIMPLE Algorithm

AortaCFD uses the PIMPLE algorithm Issa 1986, which merges the pressure-implicit splitting of PISO (Pressure-Implicit with Splitting of Operators) with the iterative relaxation of SIMPLE (Semi-Implicit Method for Pressure-Linked Equations). The outer correction loop of PIMPLE enables larger time steps than pure PISO while preserving temporal accuracy through convergence-based iteration.

4.1 Algorithm Structure

Each time step proceeds as follows:

1. Outer corrector loop (up to nOuterCorrectors iterations):
   1. Assemble and solve the momentum equation (momentum predictor).
   2. Inner corrector loop (up to nCorrectors pressure corrections):
      1. Construct the pressure equation from the semi-discretised momentum equation.
      2. Solve the pressure equation.
      3. Apply non-orthogonal corrections (nNonOrthogonalCorrectors passes).
      4. Correct the velocity field to satisfy continuity.
   3. Evaluate outer-loop residual convergence against residualControl.
   4. If residuals are below the target, exit the outer loop early.
2. Advance to the next time step.

4.2 PIMPLE Configuration by Profile

Parameter	robust	standard	precise
nOuterCorrectors	25	10	10
nCorrectors	3	2	2
nNonOrthogonalCorrectors	2	0	0
momentumPredictor	yes	yes	yes
Outer residual target (p)	\(10^{-3}\)	\(10^{-3}\)	\(10^{-3}\)
Outer residual target (U)	\(10^{-3}\)	\(10^{-4}\)	\(10^{-4}\)

4.3 Relaxation Factors

Factor	robust	standard	precise
\(p\)	0.3	0.5	0.5
\(p_{\text{Final}}\)	1.0	1.0	1.0
\(U\)	0.7	0.8	0.8
\(U_{\text{Final}}\)	1.0	1.0	1.0

Critical: Final-Iteration Relaxation

The settings \(p_{\text{Final}} = 1.0\) and \(U_{\text{Final}} = 1.0\) on the final outer corrector iteration are mandatory for correct Windkessel outlet coupling. Under-relaxation on the final corrector causes the Windkessel ODE state to drift from the applied pressure field, producing non-physical pressure waveforms and potential divergence. This requirement applies to all three profiles.

4.4 Convergence-Based Early Exit

All profiles set a deliberately high nOuterCorrectors ceiling combined with outerCorrectorResidualControl for convergence-based early exit. The solver iterates until the specified residual targets are reached, then exits the outer loop regardless of the remaining iteration budget. This strategy ensures adequate convergence during haemodynamically demanding phases (peak systole, flow reversal) while permitting rapid exit during quiescent diastolic periods where convergence is reached in two to three iterations.

Inner Solver Tolerances

The residualControl block in fvSolution sets the absolute tolerance for the linear solvers (typically GAMG for pressure, PBiCGStab for velocity). The robust profile uses \(10^{-4}\), standard uses \(10^{-6}\), and precise uses \(10^{-8}\).

---

5. Numerics Profiles Summary

AortaCFD ships three self-consistent numerics profiles. Each profile specifies temporal, spatial, and solver settings that have been tested together for stability and accuracy on patient-specific aortic meshes.

Characteristic	robust	standard	precise
Formal spatial order	1st	2nd	2nd
Formal temporal order	1st	2nd	2nd
Numerical diffusion	High	Low--moderate	Low
Stability	Maximum	High	Good (mesh-dependent)
Mesh requirement (ortho)	Any	> 50 deg	> 70 deg
Mesh requirement (skew)	Any	< 4	< 2
Intended use	Debugging, poor meshes	Production, clinical	LES, validation
Computational cost	1x (baseline)	1x	2--3x

Profile Selection Guidance

1. Start with standard for all production simulations.
2. If divergence occurs, switch to robust and investigate mesh quality.
3. Use precise only after checkMesh confirms orthogonality > 70 deg and skewness < 2.
4. For LES, precise is mandatory (LUST preserves resolved eddies).

---

6. Physics Models

6.1 Laminar (Recommended Default)

Aortic Reynolds numbers typically range from 500 to 4000, placing the flow in the transitional regime. Published patient-specific aortic CFD increasingly adopts laminar simulations for \(\text{Re} < 4000\), on the grounds that turbulence models developed for fully turbulent flows may introduce greater modelling error than the unresolved transitional fluctuations Steinman 2002; Taylor & Figueroa 2009. AortaCFD uses the laminar solver module by default.

6.2 RANS: k--omega SST

The \(k\text{--}\omega\) SST (Shear Stress Transport) model Menter 1994 is provided for cases with sustained turbulence (\(\text{Re} > 5000\)), such as severe aortic stenosis or mechanical heart valve flows. The transport equations are:

Turbulent kinetic energy:

\[ \frac{\partial k}{\partial t} + \nabla \cdot (\mathbf{u}\,k) = P_k - \beta^{*}\omega k + \nabla \cdot \bigl[(\nu + \sigma_k \nu_t)\,\nabla k\bigr] \]

Specific dissipation rate:

\[ \frac{\partial \omega}{\partial t} + \nabla \cdot (\mathbf{u}\,\omega) = \frac{\gamma}{\nu_t} P_k - \beta \omega^{2} + \nabla \cdot \bigl[(\nu + \sigma_\omega \nu_t)\,\nabla \omega\bigr] + 2(1 - F_1)\frac{\sigma_{\omega 2}}{\omega}\,\nabla k \cdot \nabla \omega \]

where \(P_k\) is the production of turbulent kinetic energy, \(F_1\) is the SST blending function, and the model coefficients (\(\beta^{*}\), \(\beta\), \(\gamma\), \(\sigma_k\), \(\sigma_\omega\)) are blended between the \(k\text{--}\omega\) (near-wall) and \(k\text{--}\varepsilon\) (free-stream) sets via \(F_1\) Menter 1994.

Caution: RANS in Predominantly Laminar Flow

The \(k\text{--}\omega\) SST model may over-predict eddy viscosity in predominantly laminar aortic flow, artificially increasing dissipation and damping physiologically relevant secondary flow structures Pope 2000. Users should verify that RANS results do not differ qualitatively from laminar solutions before adopting a turbulence model.

6.3 LES (Large Eddy Simulation)

For time-resolved turbulence investigations (jet breakdown distal to stenoses, aortic arch vortex dynamics), AortaCFD supports subgrid-scale (SGS) modelling:

• WALE (Wall-Adapting Local Eddy-viscosity) Nicoud & Ducros 1999 --- preferred for wall-bounded flows; the SGS viscosity naturally vanishes at solid walls without requiring damping functions.
• Smagorinsky --- classical constant-coefficient model; suitable for free-shear flows but requires van Driest damping near walls.

LES Mesh and Temporal Requirements

LES demands significantly higher resolution than RANS:

• Mesh orthogonality > 70 deg, skewness < 2
• Near-wall resolution: \(y^{+} < 1\) (mandatory wall-resolved LES)
• Boundary-layer expansion ratio < 1.3
• Maximum Courant number \(\leq 0.5\)
• Use the precise numerics profile exclusively

---

7. Convergence Monitoring

Convergence is assessed at two levels:

1. Linear solver convergence (inner iterations): monitored via absolute residual tolerances in fvSolution. Typical targets are \(10^{-4}\) (robust), \(10^{-6}\) (standard), or \(10^{-8}\) (precise).

2. PIMPLE outer-loop convergence: the outerCorrectorResidualControl block specifies per-field residual targets that trigger early exit from the outer correction loop.

Recommended post-simulation checks include:

• Initial residuals for \(p\) and \(\mathbf{U}\) should decrease monotonically within each time step.
• Global mass conservation error should remain below 0.01%.
• For pulsatile simulations, examine convergence behaviour separately during systolic acceleration, peak systole, and diastole.

Mesh Independence Verification

Regardless of the numerics profile selected, solution accuracy must be confirmed through a mesh independence study. The Grid Convergence Index (GCI) method Versteeg & Malalasekera 2007 on at least three mesh levels (refinement ratio \(r \geq \sqrt{2}\)) is the accepted standard for quantifying spatial discretisation uncertainty in cardiovascular CFD.

---

8. Mesh-Adaptive Numerics

The regenerate-numerics pipeline step automatically adjusts discretisation schemes based on mesh quality metrics reported by checkMesh. This is implemented in src/config/mesh_adaptive_solver.py.

python run_patient.py BPM120 --steps regenerate-numerics

The mesh-quality-aware adjustment classifies the mesh into one of five tiers based on maximum skewness, non-orthogonality, and aspect ratio:

Tier	Skewness	Non-Orthogonality	Aspect Ratio
Excellent	< 1.5	< 55 deg	< 20
Good	< 2.5	< 65 deg	< 30
Fair	< 4.0	< 70 deg	< 50
Poor	< 6.0	< 75 deg	< 100
Critical	\(\geq\) 6.0	\(\geq\) 75 deg	\(\geq\) 100

Solver settings (correctors, relaxation factors, tolerances) are then scaled relative to the selected numerics profile to ensure stability on the actual mesh. Users may override any adjusted parameter via config.json.

---

9. Parallel Execution

Domain decomposition is performed via the scotch algorithm (default) or hierarchical methods. The scotch partitioner minimises inter-processor communication by optimising the surface-area-to-volume ratio of each subdomain.

{
  "run_settings": {
    "solution_type": "parallel",
    "subdomains": 8,
    "decomposition_method": "scotch"
  }
}

Subdomain Count

The number of subdomains should be set equal to the number of available physical cores. Over-decomposition (more subdomains than cores) increases MPI communication overhead without computational benefit. For typical patient-specific aortic meshes (1--5 million cells), 4--16 subdomains provide efficient scaling.

---

References

Full bibliography on the References page.

Found an issue or have a suggestion for this page?

Open a GitHub issue