AortaCFD-X
A data-driven augmented finite-volume method (DD-AFVM) for patient-specific aortic hemodynamics. A 1D Navier–Stokes pulse wave solver is combined with Womersley analytical lifting and an SE(3)-equivariant, divergence-free residual network that learns only where analytical theory fails. The result is a 3D velocity field at a fraction of the cost of a direct CFD solve, with conservation enforced by construction.
Status: layers 1–5 implemented, 291 unit tests passing, the Alastruey 2007 Circle-of-Willis benchmark reproduced. Clinical validation on six real patients is ongoing; the source repository is private while the manuscript is in preparation.
1. Motivation
Direct 3D CFD of a patient-specific aorta (5–20 million finite volumes, one cardiac cycle at CFL ≈ 0.5) costs 100–500 CPU-hours on a research cluster. Uncertainty quantification — 10⁴ Monte-Carlo samples to bracket the pressure gradient across a coarctation — is prohibitive at that price.
A generic graph neural network surrogate is not a clean replacement either. The receptive field is small relative to the physics: an inlet waveform has to traverse the vessel tree through message passing, and the learned weights do not generalise to unseen vessel topology. A mesh independence study on such a surrogate is meaningless.
DD-AFVM rests on a simpler observation: most of the flow in a healthy aorta is already predicted, to leading order, by analytical theory. The axial profile inside a straight pipe under a prescribed pulsatile flow is the Womersley solution — a closed-form sum of Bessel functions. A 1D Navier–Stokes solver on the vessel centreline graph supplies the inlet boundary condition of that Womersley problem at every cross-section and every time step. Together they cover roughly 75 % of the FV cells with zero learned parameters.
The remaining 25 % is where analytical theory fails: the inner curve of the aortic arch (Dean secondary flow), the vicinity of bifurcations (skewing and splitting jets), the jet-and-recirculation zone downstream of a stenosis, and the first few diameters after the inlet. These are the cells the residual network corrects. The rest are untouched.
2. The 1D backbone
Each vessel segment is represented as a 1D conservation law in axial coordinate \(s\) and time \(t\):
Here \(A(s,t)\) is the cross-sectional area, \(Q(s,t)\) is the flow rate, and \(\alpha\) is the momentum-flux shape coefficient (\(\alpha=1\) for near-plug aortic flow). The friction coefficient is \(K_R = 2(\gamma + 2)\pi\mu/\rho\) with \(\gamma=9\) for a plug profile.
Closure is the Reymond–Laplace elastic tube law
which gives a wave speed \(c(A) = \sqrt{\beta \sqrt{A}/(2\rho)}\). For a healthy descending aorta, \(Eh/R_0 \approx 100\) kPa gives \(c \approx 7\) m/s at diastolic area.
Discretisation. Finite volume in space with MUSCL reconstruction (minmod limiter), Rusanov flux, and two-stage explicit Heun time stepping under CFL \(\le 0.5\). The code is in aortacfdx/onedim/pulse_solver.py.
Junctions and stenoses. At every branch or narrowing, three equations are enforced simultaneously: mass conservation, total-pressure (or static-pressure, for stenoses) continuity, and the Riemann invariant carried by each incoming characteristic. The resulting algebraic system is solved by Newton–Raphson at every time step. A stenosis is coupled through the Young–Tsai pressure-jump law (Yuhn 2022, Eq. 4), shown below.
Outlets carry a three-element Windkessel (RCR). The characteristic impedance \(R_1 = \rho c_0 / A_0\) is set to minimise wave reflection; \(R_2\) and \(C\) are calibrated by the Pietranski–Runge (PR) loop against the measured flow split and mean arterial pressure.
Validation. The Alastruey 2007 Circle-of-Willis benchmark (33 vessels, 14 bifurcations, 4 anastomoses) reproduces mean arterial pressure 79.4 mmHg, cardiac output 5.75 L/min, and bilateral cerebral symmetry within 0.3 %.
Interactive: pulse wave on an elastic tube
The widget runs the full 1D system in-browser with MUSCL–Rusanov flux and two-stage Heun time stepping. A Gaussian flow pulse is injected at the inlet; the tube dilates and contracts according to \(p(A)\), and pressure (colour) travels at the local wave speed \(c(A) = \sqrt{\beta\sqrt{A}/2\rho}\). Increasing \(Eh/R_0\) from 30 to 300 kPa raises \(\beta\) and visibly accelerates propagation. The outlet boundary toggles between a 3-element Windkessel and a non-reflecting condition.
Interactive: Young–Tsai stenosis
The pressure drop across an axisymmetric stenosis is the sum of three terms:
The viscous resistance \(R_v = \int_0^{L_s} 128\mu / (\pi D(x)^4)\,\mathrm{d}x\) is evaluated by trapezoidal quadrature on 201 points over a cosine-tapered diameter profile, exactly as in the implementation. The kinetic term scales as \(1/(1-\mathrm{SR})^4\) at high severity and dominates above SR ≈ 0.5.
3. Womersley lifting
Inside a straight segment of radius \(R\), the axial velocity profile under a flow \(Q(t) = \sum_n Q_n\,e^{in\omega t}\) is the Womersley solution
where \(\alpha_n = R\sqrt{n\omega/\nu}\) is the Womersley number of harmonic \(n\), and the shape factor
enforces the correct mean velocity. At aortic systole the fundamental \(\alpha_1 \approx 15\): the core is nearly uniform, and the velocity gradient lives in a Stokes layer of thickness \(\delta \approx R/\alpha \approx 0.8\) mm at \(R=12\) mm.
Interactive: Womersley profile
The left panel shows \(u_z(r,t)\) reconstructed from the exact formula above; the complex Bessel functions \(J_0, J_1\) are evaluated by power series. The right panel is the imposed \(Q(t) = Q_0 + Q_1 \cos(\omega t)\) with the current phase marked. At \(\alpha \to 0\) the profile collapses to Poiseuille's parabola; at \(\alpha \approx 15\) the core flattens and the velocity gradient concentrates in a Stokes layer near the wall, with peak wall shear lagging peak flow.
Accuracy of the lifting. Integrating the reconstructed 3D profile over the cross-section recovers \(Q(t)\) from the 1D solver to 0.5 %. The axis-independence test passes on four pipe orientations (the lifter is written in arc-length coordinates and rotates consistently with the local tangent \(\hat{e}_s\)).
4. Where does Womersley fail?
A finite volume cell \(x\) is flagged for neural correction if any of four geometric criteria is met:
| Criterion | Threshold | Physics it captures |
|---|---|---|
| local curvature | \(\kappa \cdot R > 0.1\) | Dean secondary flow (aortic arch) |
| junction distance | \(d_\mathrm{jct} < 3 D_\mathrm{local}\) | skewed/splitting flow at branches |
| stenosis distance | \(d_\mathrm{sten} < 10 D_\mathrm{throat}\) | jet, vena contracta, recirculation |
| entrance length | \(s_\mathrm{inlet} < \min(0.06\,\mathrm{Re}\,D,\, 10 D)\) | fully-developed assumption violated |
These thresholds are not tuned: they are conservative lower bounds from the fluid-mechanics literature. On a typical patient aorta with an arch and a mid-descending coarctation, 15–25 % of cells flip the flag. The remaining 75–85 % are predicted by Womersley alone, with no learned parameters. This is the central data-efficiency lever in the method.
Curvature is computed via the Menger formula on discrete centreline triples; distances are 3D Euclidean on the graph. The code is in aortacfdx/mask/correction_mask.py.
Interactive: where the network is needed
A schematic aorta sampled with finite-volume cells — ascending segment, arch with three supra-aortic branches, descending segment with an adjustable coarctation. Each dot is a cell, red where any of the four criteria fires. The four sliders reproduce the numbers from correction_mask.py exactly. With default thresholds and a moderate coarctation the red fraction sits between 15 and 25 %.
5. A residual network that cannot violate conservation
On the masked cells the network learns a correction \(\Delta U(x,t) = U_\mathrm{CFD}(x,t) - U_W(x,t)\). The architecture does not emit \(\Delta U\) directly. The network outputs a vector potential \(A(x,t)\), and the correction is reconstructed as
Three properties hold by construction, not as soft penalties:
- Divergence-free. \(\nabla \cdot (\nabla \times A) \equiv 0\) as a mathematical identity.
- No-slip. \(\psi_\mathrm{wall}(x) = \tanh(d_\mathrm{wall}/\ell_\mathrm{wall})\) with \(\ell_\mathrm{wall} = 0.5\) mm, so the correction vanishes smoothly at the vessel wall.
- Bounded magnitude. The scaling \(|U_W| + U_\mathrm{ref}\) with \(U_\mathrm{ref}=0.1\) m/s prevents unphysical amplification in low-flow regions.
The curl is computed on the finite-volume point cloud by a \(k\)-NN least-squares gradient operator (\(k = 20\), Tikhonov-regularised). The network layers are hand-rolled, SE(3)-equivariant, and restricted to scalar-gated linear combinations plus cross products of vector channels — sufficient for the structures present on the aortic geometry, without the overhead of the general \(e_3\) tensor-product machinery.
The loss is mean-squared error on the masked cells plus a diagnostic divergence penalty that should read as machine epsilon if the curl operator is implemented correctly:
Interactive: conservation by construction
Three panels: the vector potential \(A_z(x,y)\) on the left, the curl \(\Delta U = \nabla \times A\) in the middle (arrows), and the divergence \(\nabla \cdot (\nabla \times A)\) on the right (colour = magnitude). For every preset potential the divergence stays at the finite-difference truncation floor (≈ 10⁻¹⁴ to 10⁻¹⁶ in relative terms), independent of the correction magnitude. This is a property of the architecture, not of the loss.
6. Artefacts
| Artefact | Location | Content |
|---|---|---|
| 1D solver | aortacfdx/onedim/ |
MUSCL-Rusanov, 291 unit tests |
| Alastruey benchmark | examples/alastruey_validation.py |
33-vessel CoW, reproduces MAP and bilateral symmetry |
| Young–Tsai accuracy | outputs/experiment/pressure_gradient.json |
\(\Delta P\) error vs 3D CFD across severity |
| Womersley validation | outputs/experiment/fig_womersley_validation.pdf |
cross-section mass error < 0.5 % |
| Residual network | aortacfdx/residual/ |
22 tests including equivariance under rotation |
No trained weights exist yet. The next milestones are synthetic data generation (80 Blender geometries through AortaCFD-app) and clinical validation on six real patients. The numerics, benchmarks and infrastructure are in place; the manuscript is in preparation.
7. Parametric geometry model
The 80 training geometries are generated from a 16-parameter anatomical model of the aortic arch — two surface methods (swept tube and smooth loft), coarctation morphology from shelf to hourglass, and centreline extraction. A simpler Y-junction primitive is provided for testing junction-resolved kernels.
:material-cube-scan:{ .lg } Aortic Arch Builder → · :material-source-branch:{ .lg } Bifurcation Builder →
References
Full bibliography on the References page.
Found an issue or have a suggestion for this page?
Open a GitHub issue