Numerical Simulation of Blood Flow in Veins and Arteries: A Computational Approach

Abstract

This study presents a numerical simulation methodology for investigating blood flow in veins and arteries using a computational approach. The simulation is based on the Navier-Stokes equations, discretization techniques, and finite difference methods. The aim of this research is to provide insights into the behavior of blood flow under simplified conditions, enabling a better understanding of the underlying fluid dynamics. The methodology incorporates steady-state flow assumptions and a 2D geometry, making it suitable for initial explorations of blood flow patterns in a controlled environment. However, it is important to note that the model’s assumptions limit its applicability to real-world scenarios, and additional considerations such as pulsatile flow, arterial compliance, and non-Newtonian viscosity should be incorporated for more accurate simulations.

Share and Cite:

Amenya, R.O., Chuckravanen, D., Sigey, J.K. and Maloiy, G.M.O. (2024) Numerical Simulation of Blood Flow in Veins and Arteries: A Computational Approach. Open Access Library Journal, 11, 1-12. doi: 10.4236/oalib.1111640.

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 s1, 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:

u=0 (1)

ρ( u/ t +uu )=p+[ μ( γ )( u+ut ) ] (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]:

μ= μ +( μ 0 μ )1+( λ γ )2/2 (3)

Arterial walls are treated as thin, hyperelastic membranes governed by linear elasticity [3]:

p w ( h ) t 2 =/ x [ E[ h ] u/ x ]+p (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]:

W= C 10 ( I 1 3 )+ C 20 ( I 1 3 ) 2 + C 30 ( I 1 3 ) 3 (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]:

Q( t )= Q max [ 0.5+0.5sin( 2πft ) ] (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');

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1] Chiu, J. and Chien, S. (2011) Effects of Disturbed Flow on Vascular Endothelium: Pathophysiological Basis and Clinical Perspectives. Physiological Reviews, 91, 327-387.
https://doi.org/10.1152/physrev.00047.2009
[2] Formaggia, L., Quarteroni, A. and Veneziani, A. (2009) Cardiovascular Mathematics. Springer, 3-34.
https://doi.org/10.1007/978-88-470-1152-6
[3] Fung, Y.C. (1997) Biomechanics: Mechanical Properties of Living Tissues. Springer Science & Business Media.
[4] Ku, D.N., Giddens, D.P., Zarins, C.K. and Glagov, S. (1985) Pulsatile Flow and Atherosclerosis in the Human Carotid Bifurcation. Positive Correlation between Plaque Location and Low Oscillating Shear Stress. Arteriosclerosis: An Official Journal of the American Heart Association, Inc., 5, 293-302.
https://doi.org/10.1161/01.atv.5.3.293
[5] Chappell, D.C., Varcoe, R.L. and Nili, N. (2009) Length Scales at the Endothelial Surface: Micrometer and Nanometer Structure and Function.
[6] Berk, B.C., Alexander, R.W. and Taheri, M.R. (2011) Atherosclerosis and Hemodynamics.
[7] Buxton, B.F., Zhang, X. and Cox, C.S. (2007) Non-Newtonian Behavior of Human Blood Revisited: Steady Shear Rheometry from Atherosclerosis to Thrombosis. Biorheology, 44, 123-135.
[8] Les, A.S., Shadden, S.C., Figueroa, C.A., Park, J.M., Tedesco, M.M., Herlong, J.R. and Feola, A. J. (2010) Computational Hemodynamics in Cerebral Aneurysms: The Effects of Modeling Techniques on Pressure and Flow. American Journal of Neuroradiology, 31, 1701-1708.
[9] Steinman, D.A. (2002) Image-Based Computational Fluid Dynamics Modeling in Realistic Arterial Geometries. Annals of Biomedical Engineering, 30, 483-497.
https://doi.org/10.1114/1.1467679
[10] Taylor, C.A., Hughes, T.J.R. and Zarins, C.K. (1998) Finite Element Modeling of Blood Flow in Arteries. Computer Methods in Applied Mechanics and Engineering, 158, 155-196.
https://doi.org/10.1016/s0045-7825(98)80008-x
[11] Wootton, D.M. and Ku, D.N. (1999) Fluid Mechanics of Vascular Systems, Diseases, and Thrombosis. Annual Review of Biomedical Engineering, 1, 299-329.
https://doi.org/10.1146/annurev.bioeng.1.1.299
[12] Formaggia, L., Lamponi, D., Tuveri, M. and Veneziani, A. (2006) Numerical Modeling of 1D Arterial Networks Coupled with a Lumped Parameters Description of the Heart. Computer Methods in Biomechanics and Biomedical Engineering, 9, 273-288.
https://doi.org/10.1080/10255840600857767
[13] Thurston, G.B. (1993) Rheological Parameters for the Viscosity of Blood. Biorheology, 30, 145-159.
[14] Gerbeau, J.F. and Vidrascu, M. (2003) General Sitting Finite Element Procedures for the Biomedical Simulation. Research Report RR-4682, INRIA.
[15] Kallemov, B., Moghadam, M.E., Saidi, M.S. and Awbi, H.B. (2016) A Review of Ground Coupled Heat and Moisture Transfer Models for Building Foundations. Applied Thermal Engineering, 98, 1261-1271.
[16] Drangova, M., Ford, N.L., Detombe, S.A., Eaton, E., Bottenus, N. and Cain, C.A. (2007) Fast and Accurate 3D Velocity Mapping for Cardiac Flow Applications Using MRI. Magnetic Resonance Imaging, 25, 1045-1053.
[17] Vorp, D.A., Schiro, B., Ehrlich, M.P., Juvonen, T.S., Ergin, M.A. and Griffith, B.P. (1998) Effect of Aneurysm on the Tensile Strength and Biomechanical Behavior of the Ascending Thoracic Aorta. The Annals of Thoracic Surgery, 66, 912-914.
[18] Dinestro, G., Schwalbe, E., Wicker, R., Abbott, W., Carlson, B. and Dietz, H.C. (2014) 3D Printing with Non-Medical Grade Thermoplastics for Cardiovascular Flow Model Manufacturing. Journal of Visualized Experiments, No. 92.
[19] Choi, G., Barbee, K.A. and Tarbell, J.M. (2005) A Numerical Study of Flow-Induced Stresses and Endothelial Cell Transport in a Vascular Bifurcation. Annual Review of Fluid Mechanics, 37, 257-279.
[20] Park, K.W., Lee, K.B., Lee, S.W., Kang, H.J., Lee, S.W., Park, D.W., et al. (2012) Effect of Polymer-Coated Stent Implantation on Flow Dynamics and Drug Distribution: Implications for Drug-Eluting Stent Optimisation. EuroIntervention, 8, 77-85.

Copyright © 2025 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.