Tutorial 6: Advanced -- Mesh Convergence and LES

Goal: Perform publication-quality verification studies: mesh convergence with GCI analysis, and Large Eddy Simulation for transitional aortic flow.

Hands-On Tools

Download: generate_methods_paragraph.py (5 KB) Time: instant | Requires: Python 3 only

Generate a copy-paste methods paragraph for your manuscript:

python3 generate_methods_paragraph.py --cpd 15 --profile standard --bp 120/80 --cycles 3 --cells 1950000

Prerequisites

• Completed Tutorials 1--5 (full AortaCFD workflow including patient-specific simulation)
• Familiarity with hemodynamic metrics (TAWSS, OSI, RRT, pressure drop)
• OpenFOAM 12 installed and sourced
• Python environment with AortaCFD dependencies installed
• At least one completed patient simulation

---

Learning Objectives

By the end of this tutorial you will be able to:

1. Perform a three-level mesh convergence study and compute the Grid Convergence Index.
2. Assess whether numerical results are independent of mesh resolution.
3. Configure and understand the requirements for Large Eddy Simulation in aortic flows.
4. Generate a methods paragraph suitable for journal submission.
5. Evaluate sensitivity of hemodynamic metrics to numerical scheme selection.

---

1. Mesh Convergence Study

1.1 Why Mesh Convergence Matters

Every finite-volume simulation introduces spatial discretisation error. The magnitude of this error depends on the mesh resolution and the order of the numerical scheme. Without demonstrating that results are independent of the mesh, numerical findings cannot be considered reliable. The Grid Convergence Index (GCI) provides a systematic, quantitative framework for assessing this error Roache 1994.

Publication Requirement

Mesh convergence evidence is a minimum reporting requirement for computational haemodynamics publications Steinman & Migliavacca 2018. Reviewers will reject manuscripts that present WSS values without demonstrating mesh independence. At minimum, two mesh levels are required; three levels are preferred.

1.2 Running a Three-Level Study

AortaCFD controls mesh resolution through the cells_per_diameter (cpd) parameter. A three-level study uses a consistent refinement ratio between levels:

# Coarse mesh (cpd = 10)
python run_patient.py BPM120 \
    --config config_tutorial_coarse.json \
    --run-name gci_coarse \
    --steps case,mesh,boundary,solver,postprocess

# Medium mesh (cpd = 15)
python run_patient.py BPM120 \
    --config config_cpd15.json \
    --run-name gci_medium \
    --steps case,mesh,boundary,solver,postprocess

# Fine mesh (cpd = 22)
python run_patient.py BPM120 \
    --config config_cpd22.json \
    --run-name gci_fine \
    --steps case,mesh,boundary,solver,postprocess

Level	cells/D	Expected Cell Count	Approximate Runtime (4 cores)
Coarse	10	100--400 k	~2 min
Medium	15	500 k -- 2 M	~15 min
Fine	22	2--5 M	~60 min

Configuration Files

Create the cpd-specific configuration files by copying config_tutorial_coarse.json and changing only the mesh.cells_per_diameter value. Keep all other parameters (boundary conditions, numerics profile, simulation time) identical to isolate the effect of mesh resolution.

1.3 Richardson Extrapolation and GCI

The GCI for the fine mesh solution is defined as:

\[ \text{GCI}_{\text{fine}} = \frac{F_s \cdot |\epsilon|}{r^p - 1} \]

where:

• \(F_s = 1.25\) is the safety factor for a three-grid study,
• \(\epsilon = (f_{\text{medium}} - f_{\text{fine}}) / f_{\text{fine}}\) is the relative change in the quantity of interest,
• \(r = h_{\text{medium}} / h_{\text{fine}}\) is the grid refinement ratio, and
• \(p\) is the observed order of convergence.

The observed order of convergence is computed from:

\[ p = \frac{\ln \left| (f_{\text{coarse}} - f_{\text{medium}}) / (f_{\text{medium}} - f_{\text{fine}}) \right|}{\ln(r)} \]

For second-order numerical schemes (the standard profile), the theoretical order is \(p = 2\). The observed order may differ due to non-uniform mesh refinement, boundary layer effects, and geometric complexity.

1.4 Interpreting GCI Results

Metric	Target GCI	Notes
Pressure drop	< 5%	Converges rapidly; typically achieved at medium resolution
TAWSS (mean)	< 10%	Converges more slowly due to near-wall gradient sensitivity
TAWSS (P95/P99)	< 15%	Percentile metrics are inherently more variable
OSI	Report range	Often non-monotonic convergence; report the range across mesh levels

Non-Monotonic WSS Convergence

Wall shear stress frequently exhibits non-monotonic convergence behaviour because the near-wall velocity gradient depends on both the cell size and the boundary layer topology, which change non-uniformly with global mesh refinement. When convergence is non-monotonic, the GCI formula is not applicable. In such cases, report the range of values across mesh levels rather than a single GCI value.

1.5 When to Perform a Convergence Study

• Before any publication: mandatory for at least one representative case in the cohort.
• When changing mesh parameters significantly: new cells_per_diameter, boundary layer settings, or geometry type.
• When WSS values differ substantially between profiles: this may indicate insufficient resolution rather than a genuine scheme effect.

For large cohort studies, it is acceptable to perform a detailed three-level GCI on one representative geometry and then apply the validated mesh settings across the cohort, with spot-checks on a subset of cases.

---

2. Large Eddy Simulation

2.1 When Laminar Assumptions Break Down

The default AortaCFD configuration assumes laminar flow, which is appropriate for many aortic simulations where the Reynolds number remains below the transitional regime. However, laminar assumptions become inadequate when:

• Peak systolic Reynolds number exceeds approximately 4000.
• Complex geometry (severe coarctation, post-surgical anatomy) produces persistent turbulent structures.
• Post-stenotic jet breakdown and vortex shedding are of scientific interest.
• Validation against turbulence-resolved 4D Flow MRI is required.

In these situations, Large Eddy Simulation (LES) resolves the energy-containing eddies directly and models only the smallest subgrid scales, providing a physically faithful representation of transitional and turbulent flow features that RANS cannot capture.

2.2 AortaCFD LES Configuration

AortaCFD supports LES through the precise numerics profile with the WALE subgrid-scale model Nicoud & Ducros 1999. When LES is selected, the framework automatically configures:

• WALE subgrid model (\(C_w = 0.325\))
• Precise numerics profile (Crank-Nicolson 0.9 temporal, LUST spatial)
• Maximum Courant number 0.5
• Eddy viscosity (\(\nu_t\)) initial field
• Appropriate outlet boundary conditions for backflow

{
  "physics": {
    "model": "les",
    "simulation_type": "LES",
    "les_model": "WALE"
  },
  "numerics": {
    "profile": "precise"
  },
  "simulation_control": {
    "end_time": 7.5,
    "writeInterval": 0.01
  }
}

Why WALE

The WALE model is preferred for cardiovascular flows because it vanishes automatically at walls (correct \(y^3\) near-wall scaling), vanishes in purely laminar regions (no spurious subgrid viscosity during diastole), and handles laminar-turbulent transition without parameter tuning Nicoud & Ducros 1999. The standard Smagorinsky model is unsuitable for aortic flow due to its non-zero wall behaviour and spurious dissipation in laminar regions.

2.3 Mesh Requirements for LES

LES imposes substantially stricter mesh requirements than laminar or RANS simulations:

Parameter	Laminar	LES
cells_per_diameter	12--15	20--30
\(y^+\) at wall	N/A	< 1
Boundary layers	3--5	10--15
Expansion ratio	1.2	< 1.1
Cell aspect ratio	< 100	< 5 (near-isotropic)
Refinement levels	[1,1]--[2,3]	[n,n] uniform only

Uniform Refinement

Non-uniform refinement levels (e.g., Roache 1994; Steinman & Migliavacca 2018) create 8:1 cell volume jumps that introduce commutation errors in the LES filtering operation. For LES, always use uniform refinement levels [n,n] to maintain a smooth filter-width field.

2.4 Computational Cost

LES is 10--50 times more expensive than equivalent laminar simulations, and in some cases substantially more. The cost increase arises from three factors: finer mesh (5x more cells), smaller time step (maxCo 0.5 vs 1.0), and more cardiac cycles for statistical convergence (minimum 5, ideally 15+).

Factor	Laminar	LES
Mesh	2M cells	10M cells
Max Courant	1.0	0.5
Cycles	3	15
Total relative cost	1x	~250--750x

HPC Required

LES of patient-specific aortic geometry is not feasible on a desktop workstation. Plan for HPC resources with at least 100--200 cores and several days of wall time per case.

2.5 Best Practices

1. Statistical convergence: Simulate at least 5 cardiac cycles for converged mean velocity and WSS. Reynolds stresses and higher-order statistics require 20--50 cycles. Mean quantities converge as \(1/\sqrt{N}\).
2. Minimum mesh size: 3M cells absolute minimum; 5--10M recommended for patient-specific geometries.
3. Quality assessment: Check that \(\nu_{\text{SGS}}/\nu\) is in the range 0.1--0.5 in the bulk flow. Values exceeding 1.0 indicate insufficient mesh resolution.
4. Discard initial transient: Exclude the first 2--3 cardiac cycles from statistical analysis to remove startup effects.

---

3. SLURM Batch Submission

For HPC environments using the SLURM workload manager, AortaCFD can generate submission scripts:

python run_batch.py --cases BPM120 --slurm --cpus-per-task 32

Alternatively, use a manual SLURM script:

#!/bin/bash
#SBATCH --job-name=AortaCFD_BPM120
#SBATCH --nodes=1
#SBATCH --ntasks=32
#SBATCH --time=24:00:00
#SBATCH --mem=64G
#SBATCH --partition=compute

module load openfoam/12
source /opt/openfoam12/etc/bashrc

cd /path/to/AortaCFD
source venv/bin/activate

python run_patient.py BPM120 \
    --config config.json \
    --run-name hpc_production \
    --steps case,mesh,boundary,solver,postprocess

Pre-Mesh Locally

For large meshes, consider running the mesh generation step locally (or on a single node) and then submitting only the solver step to the HPC queue. This avoids consuming parallel resources for the serial meshing phase:

# Local: generate case and mesh
python run_patient.py BPM120 --steps case,mesh,boundary --run-name hpc_run

# HPC: run solver and post-process on existing case
python run_patient.py BPM120 --update output/BPM120/hpc_run --steps solver,postprocess

---

4. Methods Paragraph Generator

Writing a complete and accurate methods section for a computational haemodynamics manuscript is error-prone. The generate_methods_paragraph.py tool produces a copy-paste methods paragraph from simulation parameters:

python3 generate_methods_paragraph.py \
    --cpd 15 \
    --profile standard \
    --bp 120/80 \
    --cycles 3 \
    --cells 1950000

Example Output

Generated Methods Text

Simulations were performed using AortaCFD v2.1.0 with OpenFOAM 12 (Foundation). The patient-specific aortic geometry was reconstructed from contrast-enhanced CT angiography and meshed using snappyHexMesh with 15 cells per inlet diameter, yielding 1.95M hexahedral-dominant cells with 5 boundary layers (expansion ratio 1.2). The standard numerics profile was employed (backward temporal differencing, limitedLinearV spatial differencing) with adaptive time-stepping (maximum Courant number 1.0). Blood was modelled as an incompressible Newtonian fluid with density 1060 kg/m\(^3\) and dynamic viscosity 0.004 Pa s. Three-element Windkessel outlet boundary conditions were calibrated from measured blood pressure (120/80 mmHg) using Murray's law for outlet flow distribution. Three cardiac cycles were simulated; hemodynamic metrics (TAWSS, OSI, RRT) were computed from the final cycle. Percentile descriptors (P95, P99) were reported rather than maximum values to avoid mesh-dependent artefacts.

The generated paragraph includes all elements required by published reporting guidelines Steinman & Migliavacca 2018 and can be adapted to the specific journal format.

---

5. Sensitivity Analysis

5.1 Numerics Profile Comparison

A systematic comparison of the three numerics profiles on the same case reveals which hemodynamic metrics are sensitive to numerical scheme selection:

for PROFILE in robust standard precise; do
    python run_patient.py BPM120 \
        --config config_tutorial_coarse.json \
        --profile ${PROFILE} \
        --run-name sensitivity_${PROFILE}
done

5.2 Expected Sensitivity

Metric	Sensitivity to Profile	Explanation
Pressure drop	Low	Depends on global momentum balance, not near-wall gradients
TAWSS mean	Moderate	Affected by numerical diffusion smearing near-wall velocity
TAWSS P95/P99	High	Peak WSS depends strongly on resolved near-wall gradient
OSI	High	Oscillatory features damped by first-order numerical diffusion
RRT	High	Inherits sensitivity from both TAWSS and OSI

Practical Implication

Pressure drop is the most robust hemodynamic metric and can be reported with confidence even on coarse meshes with the robust profile. WSS-derived metrics (TAWSS, OSI, RRT) require at least the standard profile on an adequately resolved mesh to produce reliable values. If TAWSS P95 changes by more than 30% between robust and standard profiles, the mesh resolution is insufficient.

5.3 Minimum Reporting Checklist

Computational haemodynamics publications should include sufficient detail for reproducibility Steinman & Migliavacca 2018. The following checklist summarises the essential reporting elements:

Mesh:

• Total cell count and cells_per_diameter value
• Number of boundary layers and expansion ratio
• checkMesh quality metrics (maximum non-orthogonality, maximum skewness)
• Mesh convergence evidence (GCI or multi-level comparison)

Boundary Conditions:

• Inlet type, waveform source, and profile
• Blood pressure values used for Windkessel calibration
• Outlet model and flow distribution methodology

Simulation:

• Number of cardiac cycles and which cycle(s) used for analysis
• Time-stepping method and maximum Courant number
• Numerics profile or equivalent scheme description

Post-Processing:

• Hemodynamic metrics reported and their definitions
• Percentile descriptors and masking thresholds
• Software version

---

Exercise 1: Three-Level Mesh Convergence

Perform a GCI study on BPM120 using cells_per_diameter values of 10, 15, and 22. Record the following for each level:

Level	cells/D	Cell count	Pressure drop (mmHg)	TAWSS P95 (Pa)
Coarse	10
Medium	15
Fine	22

Compute the observed order of convergence \(p\) and the GCI for both pressure drop and TAWSS P95. Is the convergence monotonic for both quantities?

Exercise 2: Profile Sensitivity

Run BPM120 with config_tutorial_coarse.json using all three profiles (robust, standard, precise). Compare:

1. Which metric shows the largest variation across profiles?
2. Which metric is most stable?
3. At what mesh resolution would you expect the profile sensitivity to diminish?

Exercise 3: Generate a Methods Paragraph

Use generate_methods_paragraph.py with the parameters from one of your completed simulations. Compare the output against the reporting checklist in Section 5.3. Does the generated paragraph include all required elements?

Discussion

A reviewer asks: "How do you know your WSS results are mesh-independent?" Using the GCI results from Exercise 1 and the profile comparison from Exercise 2, draft a two-sentence response suitable for a reviewer rebuttal letter.

---

Summary

Topic	Key Takeaway
Mesh convergence	Three-level GCI with \(F_s = 1.25\); target < 5% for pressure, < 10% for TAWSS mean
LES	WALE model, cpd 20--30, \(y^+ < 1\), minimum 5 cycles; 250--750x laminar cost
SLURM	Pre-mesh locally, submit solver to HPC; run_batch.py --slurm for automation
Methods paragraph	generate_methods_paragraph.py produces publication-ready text
Sensitivity	Pressure drop robust to profile/mesh; WSS sensitive; always verify with convergence study

---

References

Full bibliography on the References page.

Found an issue or have a suggestion for this page?

Open a GitHub issue