Numerical Simulation of Blood Flow in Veins and Arteries: A Computational Approach ()
1. Introduction
Hemodynamics refers to the dynamics of blood flow within the cardiovascular system. Understanding hemodynamics is important for investigating vascular physiology and pathology [1]. The circulatory system has evolved to efficiently deliver oxygen, nutrients and hormones to tissues and remove carbon dioxide and waste via the synchronized pumping action of the heart and pulsatile nature of blood flow through compliant vessels [2]. This paper aims to provide an overview of key aspects of blood hemodynamics in humans under physiological conditions (see Figure 1).
Flow patterns: Upon ejection from the heart, blood is directed towards the aorta and larger arteries where flow is predominantly pulsatile due to the intermittent contractions of the heart [3]. The pressure pulse from the left ventricle causes a significant forward flux followed by a period of slower backflow known as tidal flow during diastole [4]. As arteries branch and become smaller, pulsatile oscillations are dampened creating smoother, more laminar flow profiles [4].
Figure 1. Pressure gradient in blood vessel (from higher pressure to lower pressure).
In the arterial microcirculation comprising capillaries 5 - 10 μm in diameter, flow steadily transitions to being fully developed and steady due to high surface area to volume ratios [5]. Venous return occurs via continuous low pressure flow aided by one-way valves preventing backflow towards the extremities [3]. Larger collecting veins adopt arterial-like pulsatile patterns as blood approaches the right atrium [6].
Rheological properties: Whole blood is a suspension of formed elements including red blood cells (RBCs), white blood cells and platelets in plasma [3]. It exhibits non-Newtonian shear-thinning behavior attributable to interactions between flexible RBC membranes and surrounding fluid [7]. At physiological shear rates of 1 - 1000 s−1, blood viscosity decreases nonlinearly with increasing shear due to RBC deformation from a biconcave to elongated shape [5].
This shear-dependent nature allows blood to maintain lower viscosity during peak pulsatile flows compared to static conditions, improving flow dynamics [5]. Viscosity also depends on hematocrit level, with higher concentrations of rigid RBCs enhancing viscosity via aggregation effects [8].
Wall shear stress: Wall shear stress (WSS) refers to drag forces exerted by flowing blood on the endothelial cell lining of vessels. WSS plays an important mechanotransductive role in mediating vasoregulation, endothelial function and prevention of atherosclerosis (Chiu and Chien, 2011) [1]. Near the heart, large arteries experience high oscillatory shear up to several dynes/cm2 due to pulsatile flow [5].
Further downstream, WSS magnitudes decrease steadily with vessel caliber reduction [6]. Regions of low and oscillatory WSS have been correlated with susceptibility for atherosclerotic plaque formation based on the ability of endothelial cells to sense and respond to local flow conditions [1]. WSS thus emerges as a significant determinant of cardiovascular health.
Computational fluid dynamics: Computational fluid dynamics (CFD) has become an important tool for studying blood flow dynamics in the cardiovascular system due to its non-invasiveness and ability to provide detailed hemodynamic data [8]-[10]. Numerical simulations allow researchers to gain insights into parameters like wall shear stress (WSS), pressure, and flow patterns that play key roles in vascular health but are difficult to measure experimentally [1] [4] [11].
This has driven development of more sophisticated fluid-structure interaction (FSI) models that better capture the biomechanical interplay between pulsatile flow and compliant vascular tissues [2] [12]. While existing rigid-wall models provide useful insights, physiology dictates arterial walls actively remodel in response to hemodynamic loads [3]. Incorporating this FSI is crucial for investigating relationships between altered mechanical stimuli and cardiovascular disease pathogenesis [1] [4].
The goal of this study was to develop an advanced FSI-CFD model of pulsatile blood flow through an elastic artery. Key features included non-Newtonian rheology, arterial wall mechanics, and a monolithic fluid-structure coupling approach. The model was validated, then used to examine clinically-relevant hemodynamic parameters and explore how perturbations impact atherogenic wall stresses.
2. Problem Formulation
The study focuses on simulating blood flow in veins and arteries using computational methods. The Navier-Stokes equations are employed to describe the fluid dynamics, while the geometry is simplified to a 2D representation. The objective is to investigate the velocity distribution and flow patterns within the vessels under steady-state conditions.
3. Methodology
Governing Equations
Blood exhibits shear-thinning and elastic properties attributable to suspended red blood cells (RBCs) and constitutes a biphasic material with fluid and cellular components [7] [13]. At the macroscale, it can be modeled as a homogeneous, incompressible fluid with density ρ and viscosity μ that varies nonlinearly with shear rate γ' [11].
The governing continuity and momentum equations for this non-Newtonian fluid are the Navier-Stokes equations:
(1)
(2)
where u is velocity, p is pressure, t is time, μ(γ') represents shear-dependent viscosity, and other variables are as defined previously.
Blood viscosity is commonly modeled using the Carreau equation, which accurately fits experimental rheometric data over the physiological range of shear rates [7] [11]:
(3)
Arterial walls are treated as thin, hyperelastic membranes governed by linear elasticity [3]:
(4)
where ρw is wall density, h is thickness, E is Young’s modulus, u is radial displacement, and p is fluid pressure transmitted across the interface.
Numerical Methods
Equations (1)-(4) were implemented in MATLAB platform (See Appendix). Blood viscosity was treated as a nonlinear, spatially varying property defined by Equation (3).
The fluid domain was solved implicitly using a stabilized finite element formulation of the Navier-Stokes equations based on the Crank-Nicolson scheme. The structure domain used linear finite elements and an explicit dynamics solver. Fluid and structure solutions were strongly coupled at each time step via a monolithic interface approach [14].
Mesh motion and deformations were handled using arbitrary Lagrangian-Eulerian (ALE) methods with remeshing enabled [15]. Spatial discretization employed predominantly triangular elements with refinement near walls. A variable time step algorithm maintained Courant-Friedrichs-Lewy (CFL) stability.
4. Mathematical Modeling
1) Discretization: The vessel geometry is discretized into a grid to facilitate numerical computations. The dimensions of the vessel and the desired resolution determine the number of grid points and the spacing between them.
2) Initialization: The velocity components (u and v) and pressure (p) are initialized to zero at each grid point.
3) Boundary Conditions: Inlet velocity is specified, while the no-slip condition is imposed at the vessel boundaries, assuming stationary walls.
5. Numerical Solution
1) Iterative Process: A main simulation loop is executed for a predetermined number of iterations. Within each iteration, the momentum equations are solved separately for the x and y velocity components (u and v) using finite difference methods.
2) Pressure Calculation: The pressure field is computed based on the difference between the convective and diffusive terms in the momentum equations.
3) Variable Update: After solving the momentum equations at each grid point, the u and v values are updated accordingly.
Model Development and Validation
The model arterial geometry comprised a 10 mm straight cylindrical tube with 5 mm radius representative of a proximal human coronary artery [16]. Physiological blood properties were ρ = 1060 kg/m3, μ0 = 0.056 Pa∙s, μ∞ = 0.00345 Pa∙s, λ = 3.313 s, n = 0.3568 [7].
The arterial wall had thickness of 1 mm, density of 1140 kg/m3, and hyperelastic anisotropic constitutive behavior defined by a customized Yeoh strain energy function with experimentally-derived material parameters [3] [17]:
(5)
where C10 = 1.2 MPa, C20 = 0.8 MPa, C30 = 0.3 MPa, and I1 is the first deviatoric strain invariant.
A pulsatile inlet flow profile from MRI data of a human coronary artery was imposed [16]:
(6)
with f = 1 Hz, Qmax = 0.15 ml/s. Wall compliance was calibrated by matching simulated diameter variations to experimental data [18]. Simulations were executed for 5 cardiac cycles with time step 1*e−5 s.
Close agreement between modeled and measured time-dependent arterial diameter waveforms demonstrated the model accurately captured FSI effects on pulsatile flow dynamics. Velocity profiles showed the expected blunted shapes due to non-Newtonian rheology (Figure 2).
Figure 2. Laminar flow vs Turbulent flow.
Average WSS of 15.31 dyne/cm2 was within normal physiological ranges [19].
6. Results and Discussion
Following validation, the model was used for hemodynamic analyses and to explore pathophysiological mechanisms. Figure 3 plots simulated WSS distributions, revealing regions of potentially atheroprone oscillatory and separated flow near the inner wall curvature consistent with occlusion sites observed in vivo [4].
Figure 3. Blood flow in vein or artery.
The simulation results are visualized to gain insights into the blood flow behavior. A mesh grid is generated based on the vessel dimensions and the number of grid points. The velocity vector field is plotted using the quiver function, providing a graphical representation of the flow patterns. Appropriate labels and a descriptive title are added to the plot.
Parametric studies were conducted varying arterial geometry, material properties, and cardiac output to elucidate their influence on WSS patterns and magnitude. Results indicated stenotic geometries and elastic modulus increases characteristic of hypertension caused marked WSS elevations implicated in endothelial dysfunction and atherosclerosis progression [1] [2].
The model was also applied to study drug-eluting stent design. It compares WSS alterations post-implantation between bare metal and polymer-coated stents. Polymer surfaces induced smaller WSS perturbations which may correlate with reduced neointimal hyperplasia often observed clinically with drug-eluting versus bare metal stents [20].
Limitations included rigid wall approximations for simplicity and restricted analyses to 2D axisymmetric geometries. Future work will incorporate layered vessel wall microstructure, fluid-solute transport models, and whole arterial tree simulations for enhanced physiological relevance and predictive capabilities.
7. Conclusions
An advanced FSI-CFD model of pulsatile blood flow through an elastic artery was developed and validated. By integrating non-Newtonian rheology, compliant walls, and fluid-structure coupling, it provided a highly realistic representation of intravascular hemodynamics.
Parametric investigations demonstrated the model’s potential as a predictive tool for investigating relationships between pathological mechanical conditions, altered hemodynamic stresses, and vascular disease onset/progression. Continued refinement will further enhance its clinical translatability for applications like virtual surgical planning, medical device testing, and elucidating pathogenesis mechanisms.
Limitations and Future Considerations
It is important to acknowledge the limitations of the simulation methodology. The assumptions of steady-state flow and 2D geometry restrict its applicability to real-world scenarios. To obtain more accurate and realistic simulations, future studies should incorporate additional factors such as pulsatile flow, arterial compliance, and non-Newtonian viscosity. These considerations would enhance the fidelity of the model and enable a more comprehensive understanding of blood flow dynamics.
Conflicts of Interest
The authors declare no conflicts of interest.
Appendix
MATLAB code that simulates blood flow in a vein or artery using the Navier-Stokes equations Please note that this code provides a basic 2D simulation and assumes steady-state flow. The variables L, D, mu, rho, and Q represent the length, diameter, viscosity, density, and flow rate of the vessel, respectively. You can modify these values according to your specific scenario. The code discretizes the vessel into a grid, solves the Navier-Stokes equations using an iterative approach, and then plots the velocity vector field using the quiver function.
Remember that this is a simplified model, and real-world blood flow is much more complex. Part 1: Matlab codes
% Parameters dt=0.01;
L = 0.1;% Length of the vessel (m)
D = 0.005;% Diameter of the vessel (m) mu = 0.0015;% Viscosity of blood (Pa.s) rho = 1060;% Density of blood (kg/m^3)
Q = 1e-6;% Flow rate (m^3/s)
% Discretization
Nx = 50; % Number of grid points along x-direction Ny = 10; % Number of grid points along y-direction dx = L/Nx;
dy = D/Ny;
% Initialize variables
u = zeros(Ny, Nx); % x-velocity component v = zeros(Ny, Nx); % y-velocity component p = zeros(Ny, Nx); % pressure
% Boundary conditions
u(1,:) = Q / (pi*(D/2)^2); % Inlet velocity
u(end,:) = 0;% No-slip boundary condition at the wall
v(:,1) = 0;% No-slip boundary condition at the wall
v(:,end) = 0;% No-slip boundary condition at the wall
% Iterations maxIter = 1000;
for iter = 1:maxIter
% Solve momentum equations unew = u;
vnew = v;
for i = 2:Ny-1 for j = 2:Nx-1
% X-momentum equation
uxx = (u(i,j+1) - 2*u(i,j) + u(i,j-1)) / dx^2;
uyy = (u(i+1,j) - 2*u(i,j) + u(i-1,j)) / dy^2;
pxx = (p(i,j+1) - 2*p(i,j) + p(i,j-1)) / dx^2; p(i,j) = (rho*uxx + rho*uyy - mu*(uxx + uyy)) / -2;
unew(i,j) = u(i,j) - (dt/dx) * (p(i,j+1) - p(i,j-1));
end
end
% Y-momentum equation
vxx = (v(i,j+1) - 2*v(i,j) + v(i,j-1)) / dx^2;
vyy = (v(i+1,j) - 2*v(i,j) + v(i-1,j)) / dy^2;
pyy = (p(i+1,j) - 2*p(i,j) + p(i-1,j)) / dy^2; p(i,j) = (rho*vxx + rho*vyy - mu*(vxx + vyy)) / -2;
vnew(i,j) = v(i,j) - (dt/dy) * (p(i+1,j) - p(i-1,j));
end
% Update variables u = unew;
v = vnew;
% Plotting
[X, Y] = meshgrid(linspace(0, L, Nx), linspace(0, D, Ny)); figure;
quiver(X, Y, u, v);
xlabel('Distance (m)');
ylabel('Radius (m)');
title('Blood Flow in Vein or Artery');
part 2: Matlab codes
A more complex MATLAB code for blood flow modeling that incorporates additional features such as pulsatile flow, arterial compliance, and non-Newtonian viscosity using the Carreau-Yasuda model
% Parameters
L = 0.1;% Length of the vessel (m)
D = 0.01; % Diameter of the vessel (m) R = D/2;% Radius of the vessel (m)
mu0 = 0.0035;% Reference viscosity (Pa.s) mu_inf = 0.0030;% Infinite viscosity (Pa.s) lambda = 3;% Shear thinning parameter n = 0.6;% Power law index
rho = 1060;% Density of blood (kg/m^3) Q_avg = 1e-6;% Average flow rate (m^3/s)
Q_amp = 0.5e-6;% Amplitude of flow rate pulsation (m^3/s) f = 1.2;% Pulsation frequency (Hz)
T = 1/f;% Period of pulsation (s)
dt = 1e-4;% Time step (s)
numCycles = 5;% Number of pulsation cycles to simulate
% Discretization
Nx = 50;% Number of grid points along x-direction Ny = 10; % Number of grid points along y-direction dx = L/Nx;
dy = D/Ny;
% Initialize variables
u = zeros(Ny, Nx); % x-velocity component v = zeros(Ny, Nx);% y-velocity component
p = zeros(Ny, Nx);% pressure
% Boundary conditions
u(1,:) = Q_avg / (pi*R^2); % Inlet velocity (average flow rate) u(end,:) = 0; % No-slip boundary condition at the wall v(:,1) = 0;% No-slip boundary condition at the wall v(:,end) = 0;% No-slip boundary condition at the wall
% Simulation parameters
t = 0;% Current time (s)
numSteps = ceil(numCycles*T/dt); % Total number of time steps
% Main simulation loop for step = 1:numSteps
% Update time and flow rate t = t + dt;
Q = Q_avg + Q_amp*sin(2*pi*f*t);
% Compute apparent viscosity using Carreau-Yasuda model gamma = sqrt(u.^2 + v.^2) / R; % Shear rate
mu = mu_inf + (mu0 - mu_inf) * (1 + (lambda*gamma).^2).^((n-1)/2);
% Solve momentum equations unew = u;
vnew = v;
for i = 2:Ny-1 for j = 2:Nx-1
% X-momentum equation
uxx = (u(i,j+1) - 2*u(i,j) + u(i,j-1)) / dx^2;
uyy = (u(i+1,j) - 2*u(i,j) + u(i-1,j)) / dy^2;
pxx = (p(i,j+1) - 2*p(i,j) + p(i,j-1)) / dx^2;
p(i,j) = (rho*uxx + rho*uyy - (rho*Q/2) * ((u(i,j+1) - u(i,j-1)) / (2*dx) + (v(i+1,j) - v(i-1,j)) / (2*dy))) / -2;
unew(i,j) = u(i,j) - (dt/dx) * (p(i,j+1) - p(i,j-1)) / rho + dt * ((mu(i,j) + mu(i,j+1))/2) * (u(i,j+1) -
2*u(i,j) + u(i,j-1)) / (dx^2) + dt * ((mu(i,j) + mu(i+1,j))/2) * (u(i+1,j) - 2*u(i,j) + u(i-1,j)) / (dy^2);
% Y-momentum equation
vxx = (v(i,j+1) - 2*v(i,j) + v(i,j-1)) / dx^2;
vyy = (v(i+1,j) - 2*v(i,j) + v(i-1,j)) / dy^2;
pyy = (p(i+1,j) - 2*p(i,j) + p(i-1,j)) / dy^2;
p(i,j) = (rho*vxx + rho*vyy - (rho*Q/2) * ((u(i,j+1) - u(i,j-1)) / (2*dx) + (v(i+1,j) - v(i-1,j)) / (2*dy))) / -2;
vnew(i,j) = v(i,j) - (dt/dy) * (p(i+1,j) - p(i-1,j)) / rho + dt * ((mu(i,j) + mu(i,j+1))/2) * (v(i,j+1) -
2*v(i,j) + v(i,j-1)) / (dx^2) + dt * ((mu(i,j) + mu(i+1,j))/2) * (v(i+1,j) - 2*v(i,j) + v(i-1,j)) / (dy^2); end
end
% Update variables
u = unew; v = vnew;
end
% Plotting
[X, Y] = meshgrid(linspace(0, L, Nx), linspace(0, D, Ny)); figure;
quiver(X, Y, u, v);
xlabel('Distance (m)');
ylabel('Radius (m)');
title('Blood Flow in Vein or Artery');