BloodFlowTrixi.jl
BloodFlowTrixi.jl is a Julia package that implements one-dimensional (1D) and two-dimensional (2D) blood flow models for arterial circulation. These models are derived from the Navier-Stokes equations and were developed as part of my PhD research in applied mathematics, focusing on cardiovascular pathologies such as aneurysms and stenoses.
Description
This package provides:
- 1D Blood Flow Model: This model describes blood flow along compliant arteries in a single spatial dimension. It was derived under the assumption of axisymmetric flow and accounts for arterial compliance, inertia, and frictional losses.More details about this model can be found in my corresponding publication: [Article 1D]
\[\left\{\begin{aligned} + \frac{\partial a}{\partial t} + \frac{\partial}{\partial x}(Q) &= 0 \\ + \frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}\left(\frac{Q^2}{A} + A P(a)\right) &= P(a) \frac{\partial A}{\partial x} - 2 \pi R k \frac Q {A}\\ + P(a) &= P_{ext} + \frac{Eh\sqrt{\pi}}{1-\xi^2}\frac{\sqrt{A} - \sqrt{A_0}}{A_0} \\ + R &= \sqrt{\frac{A}{\pi}} +\end{aligned}\right.\]
2D Blood Flow Model: The 2D model extends the Navier-Stokes equations under the thin-artery assumption, allowing for simulations in complex arterial geometries using curvilinear coordinates. It captures both longitudinal and angular dynamics, making it more accurate than classical 1D models while being less computationally expensive than full 3D models.
This model is described in detail in: [Article 2D]
\[\left\{\begin{aligned} + \frac{\partial a}{\partial t} + \frac{\partial}{\partial \theta}\left( \frac{Q_{R\theta}}{A} \right) + \frac{\partial}{\partial s}(Q_s) &= 0 \\ + \frac{\partial Q_{R\theta}}{\partial t} + \frac{\partial}{\partial \theta}\left(\frac{Q_{R\theta}^2}{2A^2} + A P(a)\right) + \frac{\partial}{\partial s}\left( \frac{Q_{R\theta}Q_s}{A} \right) &= P(a) \frac{\partial A}{\partial \theta} - 2 R k \frac{Q_{R\theta}}{A} + \frac{2R}{3} \mathcal{C}\sin \theta \frac{Q_s^2}{A} \\ + \frac{\partial Q_{s}}{\partial t} + \frac{\partial}{\partial \theta}\left(\frac{Q_{R\theta} Q_s}{A^2} \right) + \frac{\partial}{\partial s}\left( \frac{Q_s^2}{A} - \frac{Q_{R\theta}^2}{2A^2} + A P(a) \right) &= P(a) \frac{\partial A}{\partial s} - R k \frac{Q_s}{A} - \frac{2R}{3} \mathcal{C}\sin \theta \frac{Q_s Q_{R\theta}}{A^2} \\ + P(a) &= P_{ext} + \frac{Eh}{\sqrt{2}\left(1-\xi^2\right)}\frac{\sqrt{A} - \sqrt{A_0}}{A_0} \\ + R &= \sqrt{2A} +\end{aligned}\right.\]
Both models were designed to be used with Trixi.jl, a flexible and high-performance framework for solving systems of conservation laws using the Discontinuous Galerkin (DG) method.
Features
- 1D and 2D models for arterial blood flow.
- Derived from the Navier-Stokes equations with appropriate assumptions for compliant arteries.
- To be used with Trixi.jl for DG-based numerical simulations.
- Support for curvilinear geometries and compliant wall dynamics.
Installation
To install BloodFlowTrixi.jl, use the following commands in Julia:
julia> ]
+pkg> add Trixi
+pkg> add https://github.com/your-repo/BloodFlowTrixi.jlFuture Plans
short term
- Add second order 1D model.
- Design prim variables for 1D and 2D models.
- Add proper tests for 1D and 2D models.
- Add 3D representations of the solutions for 1D and 2D models.
- Design easy to use interfaces for users to define their own initial and boundary conditions and source terms.
long term
- Add 3D fluid-structure interaction models for complex arterial geometries.
- Design support for artery networks and simulate vascular networks using the 2D and 1D model.
- Autodiff support for 1D and 2D models for parameter optimization.
License
This package is licensed under the MIT license.
Acknowledgments
This package was developed as part of my PhD research in applied mathematics, focusing on mathematical modeling and numerical simulation of blood flow in arteries. Special thanks to the developers of Trixi.jl, whose framework was invaluable in implementing and testing these models.
BloodFlowTrixi.BloodFlowTrixiBloodFlowTrixi.BloodFlowEquations1DBloodFlowTrixi.BloodFlowEquations2DTrixi.DissipationLocalLaxFriedrichsBloodFlowTrixi.boundary_condition_outflowBloodFlowTrixi.boundary_condition_pressure_inBloodFlowTrixi.boundary_condition_slip_wallBloodFlowTrixi.flux_nonconservativeBloodFlowTrixi.frictionBloodFlowTrixi.initial_condition_simpleBloodFlowTrixi.inv_pressureBloodFlowTrixi.pressureBloodFlowTrixi.pressure_derBloodFlowTrixi.radiusBloodFlowTrixi.source_term_simpleTrixi.cons2primTrixi.fluxTrixi.initial_condition_convergence_testTrixi.max_abs_speed_naiveTrixi.prim2consTrixi.source_terms_convergence_test
BloodFlowTrixi.BloodFlowTrixi — ModulePackage BloodFlowTrixi v0.0.1This package implements 1D and 2D blood flow models for arterial circulation using Trixi.jl, enabling efficient numerical simulation and analysis.
Docs under https://yolhan83.github.io/BloodFlowTrixi.jl
BloodFlowTrixi.BloodFlowEquations1D — TypeBloodFlowEquations1D(;h,rho=1.0,xi=0.25,nu=0.04)Blood Flow equations in one space dimension. This model describes the dynamics of blood flow along a compliant artery using one-dimensional equations derived from the Navier-Stokes equations. The equations account for conservation of mass and momentum, incorporating the effect of arterial compliance and frictional losses.
The governing equations are given by
\[\left\{\begin{aligned} + \frac{\partial a}{\partial t} + \frac{\partial}{\partial x}(Q) &= 0 \\ + \frac{\partial Q}{\partial t} + \frac{\partial}{\partial x}\left(\frac{Q^2}{A} + A P(a)\right) &= P(a) \frac{\partial A}{\partial x} - 2 \pi R k \frac Q {A}\\ + P(a) &= P_{ext} + \frac{Eh\sqrt{\pi}}{1-\xi^2}\frac{\sqrt{A} - \sqrt{A_0}}{A_0} \\ + R &= \sqrt{\frac{A}{\pi}} +\end{aligned}\right.\]
BloodFlowTrixi.BloodFlowEquations2D — TypeBloodFlowEquations2D(;h,rho=1.0,xi=0.25)Defines the two-dimensional blood flow equations derived from the Navier-Stokes equations in curvilinear coordinates under the thin-artery assumption. This model describes the dynamics of blood flow along a compliant artery in two spatial dimensions (s, θ).
Parameters
h::T: Wall thickness of the artery.rho::T: Fluid density (default 1.0).xi::T: Poisson's ratio (default 0.25).nu::T: Viscosity coefficient.
The governing equations account for conservation of mass and momentum, incorporating the effects of arterial compliance, curvature, and frictional losses.
\[\left\{\begin{aligned} + \frac{\partial a}{\partial t} + \frac{\partial}{\partial \theta}\left( \frac{Q_{R\theta}}{A} \right) + \frac{\partial}{\partial s}(Q_s) &= 0 \\ + \frac{\partial Q_{R\theta}}{\partial t} + \frac{\partial}{\partial \theta}\left(\frac{Q_{R\theta}^2}{2A^2} + A P(a)\right) + \frac{\partial}{\partial s}\left( \frac{Q_{R\theta}Q_s}{A} \right) &= P(a) \frac{\partial A}{\partial \theta} - 2 R k \frac{Q_{R\theta}}{A} + \frac{2R}{3} \mathcal{C}\sin \theta \frac{Q_s^2}{A} \\ + \frac{\partial Q_{s}}{\partial t} + \frac{\partial}{\partial \theta}\left(\frac{Q_{R\theta} Q_s}{A^2} \right) + \frac{\partial}{\partial s}\left( \frac{Q_s^2}{A} - \frac{Q_{R\theta}^2}{2A^2} + A P(a) \right) &= P(a) \frac{\partial A}{\partial s} - R k \frac{Q_s}{A} - \frac{2R}{3} \mathcal{C}\sin \theta \frac{Q_s Q_{R\theta}}{A^2} \\ + P(a) &= P_{ext} + \frac{Eh}{\sqrt{2}\left(1-\xi^2\right)}\frac{\sqrt{A} - \sqrt{A_0}}{A_0} \\ + R &= \sqrt{2A} +\end{aligned}\right.\]
Trixi.DissipationLocalLaxFriedrichs — Method(dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, eq::BloodFlowEquations1D)Calculates the dissipation term using the Local Lax-Friedrichs method.
Parameters
u_ll: Left state vector.u_rr: Right state vector.orientation_or_normal_direction: Orientation or normal direction.eq: Instance ofBloodFlowEquations1D.
Returns
Dissipation vector.
BloodFlowTrixi.boundary_condition_outflow — Methodboundary_condition_outflow(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)Implements the outflow boundary condition, assuming that there is no reflection at the boundary.
Parameters
u_inner: State vector inside the domain near the boundary.orientation_or_normal: Normal orientation of the boundary.direction: Integer indicating the direction of the boundary.x: Position vector.t: Time.surface_flux_function: Function to compute flux at the boundary.eq: Instance ofBloodFlowEquations1D.
Returns
Computed boundary flux.
BloodFlowTrixi.boundary_condition_pressure_in — Methodboundary_condition_pressure_in(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)Implements a pressure inflow boundary condition where the inflow pressure varies with time.
Parameters
u_inner: State vector inside the domain near the boundary.orientation_or_normal: Normal orientation of the boundary.direction: Integer indicating the boundary direction.x: Position vector.t: Time scalar.surface_flux_function: Function to compute flux at the boundary.eq: Instance ofBloodFlowEquations1D.
Returns
Computed boundary flux with inflow pressure specified by:
\[P_{in} = \begin{cases} +2 \times 10^4 \sin^2(\pi t / 0.125) & \text{if } t < 0.125 \\ +0 & \text{otherwise} +\end{cases}\]
The corresponding inflow area A_{in} is computed using the inverse pressure relation, and the boundary state is constructed accordingly.
BloodFlowTrixi.boundary_condition_slip_wall — Methodboundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, eq::BloodFlowEquations1D)Implements a slip wall boundary condition where the normal component of velocity is reflected.
Parameters
u_inner: State vector inside the domain near the boundary.orientation_or_normal: Normal orientation of the boundary.direction: Integer indicating the direction of the boundary.x: Position vector.t: Time.surface_flux_function: Function to compute flux at the boundary.eq: Instance ofBloodFlowEquations1D.
Returns
Computed boundary flux at the slip wall.
BloodFlowTrixi.flux_nonconservative — Methodflux_nonconservative(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)Computes the non-conservative flux for the model, used for handling discontinuities in pressure.
Parameters
u_ll: Left state vector.u_rr: Right state vector.orientation::Integer: Orientation index.eq: Instance ofBloodFlowEquations1D.
Returns
Non-conservative flux vector.
BloodFlowTrixi.friction — Methodfriction(u, x, eq::BloodFlowEquations1D)Calculates the friction term for the blood flow equations, which represents viscous resistance to flow along the artery wall.
Parameters
u: State vector containing cross-sectional area and flow rate.x: Position along the artery.eq::BloodFlowEquations1D: Instance of the blood flow model.
Returns
Friction coefficient as a scalar.
BloodFlowTrixi.initial_condition_simple — Methodinitial_condition_simple(x, t, eq::BloodFlowEquations1D; R0=2.0)Generates a simple initial condition with a specified initial radius R0.
Parameters
x: Position vector.t: Time scalar.eq::BloodFlowEquations1D: Instance of the blood flow model.R0: Initial radius (default:2.0).
Returns
State vector with zero initial area perturbation, zero flow rate, constant elasticity modulus, and reference area computed as A_0 = \pi R_0^2.
This initial condition is suitable for basic tests without complex dynamics.
BloodFlowTrixi.inv_pressure — Methodinv_pressure(p, u, eq::BloodFlowEquations1D)Computes the inverse relation of pressure to cross-sectional area.
Parameters
p: Pressure.u: State vector.eq: Instance ofBloodFlowEquations1D.
Returns
Cross-sectional area corresponding to the given pressure.
BloodFlowTrixi.pressure — Methodpressure(u, eq::BloodFlowEquations1D)Computes the pressure given the state vector based on the compliance of the artery.
Parameters
u: State vector.eq: Instance ofBloodFlowEquations1D.
Returns
Pressure as a scalar.
BloodFlowTrixi.pressure_der — Methodpressure_der(u, eq::BloodFlowEquations1D)Computes the derivative of pressure with respect to cross-sectional area.
Parameters
u: State vector.eq: Instance ofBloodFlowEquations1D.
Returns
Derivative of pressure.
BloodFlowTrixi.radius — Methodradius(u, eq::BloodFlowEquations1D)Computes the radius of the artery based on the cross-sectional area.
Parameters
u: State vector.eq: Instance ofBloodFlowEquations1D.
Returns
Radius as a scalar.
BloodFlowTrixi.source_term_simple — Methodsource_term_simple(u, x, t, eq::BloodFlowEquations1D)Computes a simple source term for the blood flow model, focusing on frictional effects.
Parameters
u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.x: Position vector.t: Time scalar.eq::BloodFlowEquations1D: Instance of the blood flow model.
Returns
Source terms vector where:
s_1 = 0(no source for area perturbation).s_2represents the friction term given bys_2 = \frac{2 \pi k Q}{R A}.
Friction coefficient k is computed using the friction function, and the radius R is obtained using the radius function.
Trixi.cons2prim — MethodTrixi.cons2prim(u, eq::BloodFlowEquations1D)Converts the conserved variables to primitive variables.
Parameters
u: State vector.eq: Instance ofBloodFlowEquations1D.
Returns
Primitive variable vector.
Trixi.flux — MethodTrixi.flux(u, orientation::Integer, eq::BloodFlowEquations1D)Computes the flux vector for the conservation laws of the blood flow model.
Parameters
u: State vector.orientation::Integer: Orientation index for flux computation.eq: Instance ofBloodFlowEquations1D.
Returns
Flux vector as an SVector.
Trixi.initial_condition_convergence_test — Methodinitial_condition_convergence_test(x, t, eq::BloodFlowEquations1D)Generates a smooth initial condition for convergence tests of the blood flow equations.
Parameters
x: Position vector.t: Time scalar.eq::BloodFlowEquations1D: Instance of the blood flow model.
Returns
Initial condition state vector with zero initial area perturbation, sinusoidal flow rate, a constant elasticity modulus, and reference area.
Details
The returned initial condition has:
- Zero perturbation in area (
a = 0). - A sinusoidal flow rate given by
Q = sin(\pi x t). - A constant elasticity modulus
E. - A reference cross-sectional area
A_0 = \pi R_0^2forR_0 = 1.
This initial condition can be used to verify the accuracy and stability of numerical solvers.
Trixi.max_abs_speed_naive — MethodTrixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, eq::BloodFlowEquations1D)Calculates the maximum absolute speed for wave propagation in the blood flow model using a naive approach.
Parameters
u_ll: Left state vector.u_rr: Right state vector.orientation::Integer: Orientation index.eq: Instance ofBloodFlowEquations1D.
Returns
Maximum absolute speed.
Trixi.prim2cons — MethodTrixi.prim2cons(u, eq::BloodFlowEquations1D)Converts the primitive variables to conserved variables.
Parameters
u: Primitive variable vector.eq: Instance ofBloodFlowEquations1D.
Returns
Conserved variable vector.
Trixi.source_terms_convergence_test — Methodsource_terms_convergence_test(u, x, t, eq::BloodFlowEquations1D)Computes the source terms for convergence tests of the blood flow equations.
Parameters
u: State vector containing area perturbation, flow rate, elasticity modulus, and reference area.x: Position vector.t: Time scalar.eq::BloodFlowEquations1D: Instance of the blood flow model.
Returns
Source terms vector.
Details
The source terms are derived based on the smooth initial condition and friction effects:
s_1represents the source term for area perturbation and is given bys_1 = \pi t \cos(\pi x t).s_2represents the source term for the flow rate and includes contributions from spatial and temporal variations as well as friction effects.
The radius R is computed using the radius function, and the friction coefficient k is obtained using the friction function.
This function is useful for evaluating the correctness of source term handling in numerical solvers.
