**Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries**

Renat A. Sultanov and Dennis Guster *Department of Information Systems and BCRL, St. Cloud State University, St. Cloud, MN USA* 

#### **1. Introduction**

474 Fluid Dynamics, Computational Modeling and Applications

2. Transition from single-file to multi-file flow as a function of hematocrit in capillaries of various diameters is studied. RBCs are single-file in narrow tube and at low hematocrit,

4. Particle simulations were applied to the capillary vessel flow at finger tip. In case of bent tube, the RBCs is initially parachute shape at straight tube and then deforms to asymmetric parachute shape. In case of bent and twisted tube, initially RBC is parachute shape and then deforms to asymmetric shape (between flat and parachute

Wada, Kobayashi, Takahashi and Karino (2000). A numerical simulation of the deformation

Tanaka N., Takano T. and Masuzawa T. (2004). 3-dimensional micro-simulation of blood flow with SPH method, *Japanese Fluid Engineering Conference*, JSME, No. 712, Japan, 2004 Gaehtgens P., Duhrssen C. and Albrecht KH. (1980). Motions, Deformation and Interaction

Koshizuka S., *Computational Fluid Dynamics* (1997). in Japanese, Baihuukan, Japan, 1997 Tsubota K., Wada S. and Yamaguchi T. (2006). Particle method for computer simulation of

Nagayama K. and Tanaka K. (2004). Particle Simulations of Two Phase Blood Flow with

Nagayama K. and Tanaka K. (2005). Particle Simulations of three dimensional blood flow

Nagayama K. (2006), Particle Simulations of Blood flow in Vein with Many RBCs,

Nagayama K. and Honda K. (2008b), Particle Simulations of Blood Flow in Bent and

of an erythrocyte, *Japanese Mechanical Engineering Congress*, MECJ-05, No.1226

of Blood Cells and Plasma During Flow Through Narrow Capillary Tubes, *Blood* 

red blood cell motion in blood flow, *Computer Methods and Programs in Biomedicine*,

Red Blood Cell, *Japanese Fluid Engineering Conference*, JSME, No.G808, Japan, 2004,

with a blood cell, *Proceedings of 2005 Annual Meeting*, Japan Society of Fluid

*International Proceedings by Medimond from World Congress of Biomechanics*, pp. 557- 562 Munich, Germany, 2006, Volume ISBN 88-7587-270-8, CD ISBN 88-7587-271-6 Nagayama K. and Honda K. (2008a), Particle Simulations of the deformation of red blood

cells in a capillary vessel. *Proceedings of the 12th Asian Congress of Fluid Mechanics*,

Twisted Capillary Vessel with Red Blood Cells, *Proceedings of the TFEC*, Sapporo,

while they are multi-file as the tube diameter increases or hematocrit increases. 3. RBC shape change in time was studied for a case of ID= 7.37μm Ht = 0.49 in details. At first, RBCs were flowing in line like a parachute. But after that the shape fluctuated gradually, and it became stable in this case in the zipper state finally. A stable state is expected that it also changes by the pipe diameter, the number of erythrocytes and the

physical properties of the RBCs.

pp.287-288, Japan, 2000

*cells* Vol.6, 799-812, 1980

Vol. 83, pp. 139-146, 2006.

Mechanics AM05-17-007, Japan, 2005

Daejeon, Korea, 18-21 August 2008

Japan, 14-16 September 2008

ISSN 1348-2882

Monaghan J. (1992). *Annu. Rev. Astrophys.*, No.30, pp.543-574, 1992

**6. References** 

shape) due to the twist and bend flow.

Cardiovascular diseases, such as ischemic heart disease, myocardial infarction, and stroke are leading causes of death in Western countries. All of these vascular diseases share a common element: atherosclerosis. They also share a common final event: the failure or destruction of the vascular wall structure, Dhein et al. (2005); Waite (2005).

Atherosclerosis reduces arterial lumen size through plaque formation and arterial wall thickening. It occurs at specific arterial sites. This phenomenon is related to hemodynamics and to wall shear stress (WSS) distribution, Fung (1993). From the physical point of view WSS is the tangential drag force produced by moving blood, i.e. it is a mathematical function of the velocity gradient of blood near the endothelial surface. A general description of WSS is presented in Landau & Lifshitz (1959). Arterial wall remodeling is regulated by WSS, Grotberg & Jensen (2004), for example, in response to high shear stress arteries enlarge. From the bio-mechanical point of view one can conclude, that the atherosclerotic plaques localize preferentially in the regions of low shear stresses, but not in regions of higher shear stresses. Furthermore, decreased shear stress induces intimal thickening in vessels which have adapted to high flow.

Final vascular events that induce fatal outcomes, such as acute coronary syndrome, are triggered by the sudden mechanical disruption of an arterial wall. Thus, we can conclude, that the final consequences of tragic fatal vascular diseases are strongly connected to mechanical events that occur within the vascular wall, and these in turn are likely to be heavily influenced by alterations in blood flow and the characteristics of the blood itself.

Currently researchers in the field of biomechanics and biomedicine conduct laboratory investigations of human blood flow in different shape and size tubes, which are designed to be approximate models of human vessels and arteries, see for example Huo & Kassab (2006); Taylor & Draney (2004). Some researchers also carry out intensive computer simulations of these bio-mechanical systems, see for example Chen & Lu (2004; 2006); Cho & Kensey (1991); Duraiswamy et al. (2007); Johnston et al. (2004); Morris et al. (2004; 2005); Mukundakrishnan et al. (2008); Peskin (1977); Sultanov et al., 2008 (a;b); Sultanov & Guster (2009).

Also, there have been laboratory experiments in which specific stents are incorporated in such artificial vessels (tubes). Stent implantation has been used to open diseased coronary blood vessels, allowing improved perfusion of the cardiac muscle. Used in combination with drug therapy, vascular repair and dilation techniques (angioplasty) the implantation of metallic

**2.1 Equations**

*ρ ∂vi <sup>∂</sup><sup>t</sup>* <sup>+</sup> *vk*

simpler in form:

still open to question.

*∂vr <sup>∂</sup><sup>t</sup>* <sup>+</sup> *vr*

*∂v<sup>θ</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup> *vr* *∂vr ∂r* + *vθ r ∂vr ∂θ* <sup>+</sup> *vz*

> *∂vz ∂r* + *vθ r ∂vz ∂θ* <sup>+</sup> *vz*

*∂v<sup>θ</sup> ∂r* + *vθ r ∂v<sup>θ</sup> ∂θ* <sup>+</sup> *vz*

*∂vz <sup>∂</sup><sup>t</sup>* <sup>+</sup> *vr*

way, Landau & Lifshitz (1959):

*∂vi ∂xk*

*ρ ∂v* <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>p</sup> ∂xi* + *∂ ∂xk*

the Eq. (1) becomes the well known Navier-Stokes equation:

*<sup>∂</sup><sup>t</sup>* + (*<sup>v</sup>*∇)*<sup>v</sup>*

*<sup>∂</sup><sup>t</sup>* + (*<sup>v</sup>*∇)*<sup>v</sup>*

 *∂v*

The general form of the dynamics equation for viscous fluid can be written in the following

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 477

here *vi* and *xk* are velocity and coordinates of the fluid, *ρ* is the density of fluid, *η* and *ζ* are dynamical characteristics of the fluid, i.e. coefficients of viscosity. Because in a general case the pressure *p*, the temperature *T* and therefore the viscosity coefficients *η* and *ζ* are not constants in a flowing fluid, one cannot take them out of the partial differentials in the Eq. (1). However, in a very wide range of applications it is a good approximation to consider the variation of these coefficients to be negligible in the fluid, that is *η* = *const* and *ζ* = *const*. In these cases

Further, if the fluid is considered as incompressible: div *v* = 0, then the NS equation becomes

This equation is, probably, one of the most applicable mathematical results related to modeling and real time simulation of various physical systems, such as: air flow in aerodynamics, blood flow in medical applications and even cash flow in financial problems. The fundamental NS equation is nonlinear, diverse, rich, and, as we mentioned above, has strong practical applications in science and in various fields of modern technologies, for example, microand nano-fluidics. However, because of the considered atomistic level in the field of novel nano-fluidics problems a direct application of the NS equation to these systems might be problematic, especially at low temperatures of these systems. Further, from the mathematical point of view the NS equation is a very complicated nonlinear partial differential equation, which still remains as a "stumbling block" for mathematicians. The Clay Mathematics Institute (Boston, Massachusetts, USA) announced *Seven Millennium Problems* with a prize of US\$ 1,000,000.00 for each (http://www.claymath.org/index.php). One of these problems is related to the existence and smoothness of the NS equation. It is hard to believe that despite many successful practical applications of the NS equation its fundamental mathematical property is

Nevertheless, because this work deals with different 3D geometries it would be useful to represent the NS equation in different forms. For example, in the case of cylindrical symmetry

one can apply cylindrical coordinates (*r*, *θ*, *z*) and the NS equation could be written:

*∂vz <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂z*

*∂vr <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>v</sup>*<sup>2</sup> *θ <sup>r</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂r*

*∂v<sup>θ</sup> ∂z* + *vrv<sup>θ</sup> <sup>r</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ* 1 *r ∂p*

grad*<sup>p</sup>* <sup>+</sup> *<sup>η</sup>*

*ρ*

+ *Fr* + *ν*

*∂θ* <sup>+</sup> *<sup>F</sup><sup>θ</sup>* <sup>+</sup> *<sup>ν</sup>*

<sup>∇</sup>2*vr* <sup>−</sup> *vr*

<sup>∇</sup>2*v<sup>θ</sup>* <sup>+</sup>

*<sup>r</sup>*<sup>2</sup> <sup>−</sup> <sup>2</sup> *r*2 *∂v<sup>θ</sup> ∂θ*

> 2 *r*2 *∂vr ∂θ* <sup>−</sup> *<sup>v</sup><sup>θ</sup> r*2 , (5)

<sup>+</sup> *Fz* <sup>+</sup> *<sup>ν</sup>*∇2*vz*, (6)

, (4)

 <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ* − 2 3 *δik ∂vl ∂xl* <sup>+</sup>

= −grad*p* + *η*�*v* + (*ζ* + *η*/3)grad div*v*. (2)

*∂ ∂xi ζ ∂vl ∂xi* 

�*v*. (3)

, (1)

 *η ∂vi ∂xk* + *∂vk ∂xi*

stents has created a multibillion dollar industry. Stents are commonly used in many different blood vessels, but the primary site of deployment is in diseased coronary arteries.

Stents represent a very special case in the modeling research problems mentioned above, Frank et al. (2002). Taking into account that stents have a very small size and rather complicated structure and shape, this situation makes it difficult to obtain precise measurements. Therefore high quality and precise computer simulations of blood flow through vessels with implanted stents would be most useful, Frank et al. (2002). Work of this type is already underway and we would like to mention several pertinent studies, Benard et. al. (2006); Banerjee et. al. (2007); Faik et al. (2007); Seo et al. (2005).

Nevertheless, there are still many problems in obtaining precise realistic geometries for the required vessels. Human arteries, especially the aorta, have complicated spatial-geometric and characteristic configurations. For example, the aortic arch centerline does not lie on a plane and there are major branches at the top of the arch feeding the carotid arterial circulation. One of the main problems in the field of bio-medical blood flow simulation is to obtain precise geometrical-mathematical representations of different vessels. This information in turn needs to be included in the simulation programs.

However, it seems logical that a first step in these investigations would be to apply simpler 3D-geometry forms and models, but at the same time to take into account the precise physical effects of blood movement such as the non-Newtonian characteristics of human blood, realistic pulsatile flow, and possible turbulent effects. Because of the applied pulsatile flow in our simulations turbulence may be significant to the final results of this study.

Therefore, in the current work we carried out real-time full-dimensional computer simulations of a realistic pulsatile human blood flow in actual size vessels, vessels with a bifurcation, and in a model of the aortic arch. We take into account different physical effects, such as turbulence and the non-Newtonian nature of human blood. The next section presents the mathematical methodology and the physical model used in this work. The general purpose commercial computational fluid dynamics program FLOW3D is used for its basic functionality, but we supplemented its capability by adding our routines to obtain the results presented in this work.

Sec. 3 presents results for three vessels of different geometries. The CGS unit system is used in all simulations, as well as for presentation of the results. Conclusions and discussion comparing our results to well respected previous studies are included in Sec. 4.

#### **2. Physical models and mathematical methods**

As we mentioned above, we undertook pulsatile human blood flow simulation experiments using different size and shape human vessel arteries. For each spatial configuration one needs to provide a specific approach for the numerical solution to the complicated second order partial differential equations of the fluid dynamics. These equations are also named as the Navier-Stokes (NS) equations.

For simple cylindrical vessels we used the cylindrical coordinate system: *r* = (*r*, *θ*, *Z*). However, for the aortic arch or bifurcated vessels, where there is no cylindrical symmetry, we applied the Cartezian coordinate system: *r* = (*x*, *y*, *z*). In the cases of the aortic arch and bifurcated vessels we used up to five blocks of matched Cartezian coordinate subsystems. Below we represent the NS equations in a general form, because, for each of the special cases, considered in this work and the chosen coordinate system, the partial differential equations of the fluid dynamics may look different. At the same time we understand that the general differential operator form of these equations is unique.

#### **2.1 Equations**

2 Will-be-set-by-IN-TECH

stents has created a multibillion dollar industry. Stents are commonly used in many different

Stents represent a very special case in the modeling research problems mentioned above, Frank et al. (2002). Taking into account that stents have a very small size and rather complicated structure and shape, this situation makes it difficult to obtain precise measurements. Therefore high quality and precise computer simulations of blood flow through vessels with implanted stents would be most useful, Frank et al. (2002). Work of this type is already underway and we would like to mention several pertinent studies, Benard

Nevertheless, there are still many problems in obtaining precise realistic geometries for the required vessels. Human arteries, especially the aorta, have complicated spatial-geometric and characteristic configurations. For example, the aortic arch centerline does not lie on a plane and there are major branches at the top of the arch feeding the carotid arterial circulation. One of the main problems in the field of bio-medical blood flow simulation is to obtain precise geometrical-mathematical representations of different vessels. This information in turn needs

However, it seems logical that a first step in these investigations would be to apply simpler 3D-geometry forms and models, but at the same time to take into account the precise physical effects of blood movement such as the non-Newtonian characteristics of human blood, realistic pulsatile flow, and possible turbulent effects. Because of the applied pulsatile flow in our simulations turbulence may be significant to the final results of this study. Therefore, in the current work we carried out real-time full-dimensional computer simulations of a realistic pulsatile human blood flow in actual size vessels, vessels with a bifurcation, and in a model of the aortic arch. We take into account different physical effects, such as turbulence and the non-Newtonian nature of human blood. The next section presents the mathematical methodology and the physical model used in this work. The general purpose commercial computational fluid dynamics program FLOW3D is used for its basic functionality, but we supplemented its capability by adding our routines to obtain the results presented in this

Sec. 3 presents results for three vessels of different geometries. The CGS unit system is used in all simulations, as well as for presentation of the results. Conclusions and discussion

As we mentioned above, we undertook pulsatile human blood flow simulation experiments using different size and shape human vessel arteries. For each spatial configuration one needs to provide a specific approach for the numerical solution to the complicated second order partial differential equations of the fluid dynamics. These equations are also named as the

For simple cylindrical vessels we used the cylindrical coordinate system: *r* = (*r*, *θ*, *Z*). However, for the aortic arch or bifurcated vessels, where there is no cylindrical symmetry, we applied the Cartezian coordinate system: *r* = (*x*, *y*, *z*). In the cases of the aortic arch and bifurcated vessels we used up to five blocks of matched Cartezian coordinate subsystems. Below we represent the NS equations in a general form, because, for each of the special cases, considered in this work and the chosen coordinate system, the partial differential equations of the fluid dynamics may look different. At the same time we understand that the general

comparing our results to well respected previous studies are included in Sec. 4.

**2. Physical models and mathematical methods**

differential operator form of these equations is unique.

Navier-Stokes (NS) equations.

blood vessels, but the primary site of deployment is in diseased coronary arteries.

et. al. (2006); Banerjee et. al. (2007); Faik et al. (2007); Seo et al. (2005).

to be included in the simulation programs.

work.

The general form of the dynamics equation for viscous fluid can be written in the following way, Landau & Lifshitz (1959):

$$\rho \left( \frac{\partial v\_{\mathrm{i}}}{\partial t} + v\_{\mathrm{k}} \frac{\partial v\_{\mathrm{i}}}{\partial \mathbf{x}\_{\mathrm{k}}} \right) = -\frac{\partial p}{\partial \mathbf{x}\_{\mathrm{i}}} + \frac{\partial}{\partial \mathbf{x}\_{\mathrm{k}}} \left\{ \eta \left( \frac{\partial v\_{\mathrm{i}}}{\partial \mathbf{x}\_{\mathrm{k}}} + \frac{\partial v\_{\mathrm{k}}}{\partial \mathbf{x}\_{\mathrm{i}}} - \frac{2}{3} \delta\_{\mathrm{ik}} \frac{\partial v\_{\mathrm{l}}}{\partial \mathbf{x}\_{\mathrm{l}}} \right) \right\} + \frac{\partial}{\partial \mathbf{x}\_{\mathrm{i}}} \left( \zeta \frac{\partial v\_{\mathrm{l}}}{\partial \mathbf{x}\_{\mathrm{i}}} \right), \tag{1}$$

here *vi* and *xk* are velocity and coordinates of the fluid, *ρ* is the density of fluid, *η* and *ζ* are dynamical characteristics of the fluid, i.e. coefficients of viscosity. Because in a general case the pressure *p*, the temperature *T* and therefore the viscosity coefficients *η* and *ζ* are not constants in a flowing fluid, one cannot take them out of the partial differentials in the Eq. (1). However, in a very wide range of applications it is a good approximation to consider the variation of these coefficients to be negligible in the fluid, that is *η* = *const* and *ζ* = *const*. In these cases the Eq. (1) becomes the well known Navier-Stokes equation:

$$\rho \left[ \frac{\partial \vec{v}}{\partial t} + (\vec{v} \nabla) \vec{v} \right] = -\text{grad}\, p + \eta \triangle \vec{v} + (\zeta + \eta/3) \text{grad } \text{div}\vec{v}. \tag{2}$$

Further, if the fluid is considered as incompressible: div *v* = 0, then the NS equation becomes simpler in form:

$$
\left[\frac{\partial \vec{v}}{\partial t} + (\vec{v} \nabla) \vec{v}\right] = -\frac{1}{\rho} \text{grad}\rho + \frac{\eta}{\rho} \triangle \vec{v}.\tag{3}
$$

This equation is, probably, one of the most applicable mathematical results related to modeling and real time simulation of various physical systems, such as: air flow in aerodynamics, blood flow in medical applications and even cash flow in financial problems. The fundamental NS equation is nonlinear, diverse, rich, and, as we mentioned above, has strong practical applications in science and in various fields of modern technologies, for example, microand nano-fluidics. However, because of the considered atomistic level in the field of novel nano-fluidics problems a direct application of the NS equation to these systems might be problematic, especially at low temperatures of these systems. Further, from the mathematical point of view the NS equation is a very complicated nonlinear partial differential equation, which still remains as a "stumbling block" for mathematicians. The Clay Mathematics Institute (Boston, Massachusetts, USA) announced *Seven Millennium Problems* with a prize of US\$ 1,000,000.00 for each (http://www.claymath.org/index.php). One of these problems is related to the existence and smoothness of the NS equation. It is hard to believe that despite many successful practical applications of the NS equation its fundamental mathematical property is still open to question.

Nevertheless, because this work deals with different 3D geometries it would be useful to represent the NS equation in different forms. For example, in the case of cylindrical symmetry one can apply cylindrical coordinates (*r*, *θ*, *z*) and the NS equation could be written:

$$\frac{\partial v\_r}{\partial t} + v\_r \frac{\partial v\_r}{\partial r} + \frac{v\_\theta}{r} \frac{\partial v\_r}{\partial \theta} + v\_z \frac{\partial v\_r}{\partial z} - \frac{v\_\theta^2}{r} = -\frac{1}{\rho} \frac{\partial p}{\partial r} + F\_r + \nu \left(\nabla^2 v\_r - \frac{v\_r}{r^2} - \frac{2}{r^2} \frac{\partial v\_\theta}{\partial \theta}\right), \tag{4}$$

$$\frac{\partial v\_{\theta}}{\partial t} + v\_{r} \frac{\partial v\_{\theta}}{\partial r} + \frac{v\_{\theta}}{r} \frac{\partial v\_{\theta}}{\partial \theta} + v\_{z} \frac{\partial v\_{\theta}}{\partial z} + \frac{v\_{r} v\_{\theta}}{r} = -\frac{1}{\rho} \frac{1}{r} \frac{\partial p}{\partial \theta} + F\_{\theta} + \nu \left(\nabla^{2} v\_{\theta} + \frac{2}{r^{2}} \frac{\partial v\_{r}}{\partial \theta} - \frac{v\_{\theta}}{r^{2}}\right), \tag{5}$$

$$\frac{\partial v\_z}{\partial t} + v\_r \frac{\partial v\_z}{\partial r} + \frac{v\_\theta}{r} \frac{\partial v\_z}{\partial \theta} + v\_z \frac{\partial v\_z}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + F\_z + \nu \nabla^2 v\_z \tag{6}$$

together with the continuity equation:

$$\frac{1}{r}\frac{\partial}{\partial r}(rv\_r) + \frac{1}{r}\frac{\partial v\_\theta}{\partial \theta} + \frac{\partial v\_z}{\partial z} = 0,\tag{7}$$

where *ν* = *η*/*ρ* is the kinematic viscosity Landau & Lifshitz (1959).

However, in the general case when there is no symmetry it is useful to apply the well known Cartesian coordinates *x*, *y*, and *z*: accordingly the equations of motion for the fluid velocity components (*u*, *v*, *w*) are:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial x} + \mathbf{X} + \nu \nabla^2 u,\tag{8}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial y} + Y + \nu \nabla^2 v,\tag{9}$$

$$\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} = -\frac{1}{\rho} \frac{\partial p}{\partial z} + Z + \nu \nabla^2 w. \tag{10}$$

In this case the continuity equation has the form:

$$
\frac{\partial \mu}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0. \tag{11}
$$

0 0.2 0.4 0.6 0.8 1 Pulse Time, t [sec]

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 479

Fig. 1. The inflow time-dependent waveform used in these simulations. The function was

accelerations, (*bx*, *by*, *bz*) are the flow losses in porous media or across porous baffle plates, and the final term accounts for the injection of mass at a source represented by a geometric component. Mathematical expressions for the viscous accelerations (*fx*, *fy*, *fz*) are presented

The term *Uw* = (*uw*, *vw*, *ww*) in equations (12-14) is the velocity of the source component, which will generally be non-zero for a mass source of a General Moving Object (GMO) FLOW3D (2007). The term *Us* = (*us*, *vs*, *ws*) is the velocity of the fluid at the surface of the

where *dQ* is the mass flow rate, *ρ<sup>s</sup>* fluid source density, *dA* the area of the source surface in the cell and *n* the outward normal to the surface. The source is of the stagnation pressure type when in equations (12-14) *δ* = 0.0. Next, *δ* = 1.0 corresponds to the source of the static

As we mentioned earlier, in all simulations we considered the blood flow as a pulsatile flow. The final result for the inflow waveform has been taken from Fig. 3 of the work Papaharilaou et al. (2007). The pulse was applied for 5.5 cycle times in our work. The waveform is shown in Fig. 1. These velocity values are used as time-dependent inflow initial boundary conditions.

Next, the general mass continuity equation, which is solved within the FLOW3D program has

*∂z*

*Rdi f* is a turbulent diffusion term, and *Rsor* is a mass source. The turbulent diffusion term is

(*vpAyR∂ρ ∂y*

(*ρwAz*) + *ξ*

) + *<sup>∂</sup> ∂z* (*vpAz*

*ρuAx*

*∂ρ ∂z* ) + *ξ*

*<sup>x</sup>* <sup>=</sup> *Rdi f* <sup>+</sup> *Rsor*, (16)

*<sup>x</sup>* , (17)

*ρvpAx*

(*ρvAy*) + *<sup>∂</sup>*

*d*(*Qn*)

*dA* (15)

source relative to the source itself. It is computed in each control volume as

These numbers are included directly in the FLOW3D program.

*∂y*

) + *<sup>R</sup> <sup>∂</sup> ∂y*

(*ρuAx*) + *<sup>R</sup> <sup>∂</sup>*

*∂ρ ∂x*

(*vpAx*

*<sup>U</sup> <sup>s</sup>* <sup>=</sup> <sup>1</sup> *ρs*

in the Appendix.

pressure type.

the following *general* form:

*∂ ∂x*

*Rdi f* <sup>=</sup> *<sup>∂</sup> ∂x*

*Vf ∂ρ <sup>∂</sup><sup>t</sup>* <sup>+</sup>

taken from Fig. 3 of the work Papaharilaou et al. (2007).

Velocity, w [cm/sec]

In the specific case of the FLOW3D program the equations of motion for the fluid velocity components (*u*, *v*, *w*) with special additional terms included in the program are written:

$$\begin{split} \frac{\partial u}{\partial t} + \frac{1}{V\_F} (u A\_X \frac{\partial u}{\partial x} + v A\_Y R \frac{\partial u}{\partial y} + w A\_z \frac{\partial u}{\partial z}) - \xi \frac{A\_Y v^2}{x V\_f} &= -\frac{1}{\rho} \frac{\partial p}{\partial x} + G\_\mathbf{x} + f\_\mathbf{x} - b\_\mathbf{x} - \frac{\partial P\_f}{\partial z} \\ \frac{R\_{\rm sor}}{\rho V\_f} (u - u\_\mathbf{w} - \delta \cdot u\_\mathbf{s}) &\quad \text{(12)} \end{split} \tag{12}$$

$$\begin{split} \frac{\partial v}{\partial t} + \frac{1}{V\_F} (\mu A\_x \frac{\partial v}{\partial x} + v A\_y R \frac{\partial v}{\partial y} + w A\_z \frac{\partial v}{\partial z}) + \xi \frac{A\_y \mu v}{x V\_f} &= -\frac{R}{\rho} \frac{\partial p}{\partial y} + G\_y + f\_y - b\_y - \\ \frac{R\_{\text{sor}}}{\rho V\_f} (v - v\_w - \delta \cdot v\_s) &\quad \text{(13)} \end{split} \tag{13}$$

$$\begin{split} \frac{\partial w}{\partial t} + \frac{1}{V\_F} (u A\_X \frac{\partial w}{\partial x} + v A\_y R \frac{\partial w}{\partial y} + w A\_z \frac{\partial w}{\partial z}) &= -\frac{1}{\rho} \frac{\partial p}{\partial z} + G\_z + f\_z - b\_z - \\ &\frac{R\_{\text{sor}}}{\rho V\_f} (w - w\_w - \delta \cdot w\_s). \end{split} \tag{14}$$

Here, (*u*, *v*, *w*) are the velocity components in coordinate directions (*x*, *y*, *z*) respectively. For example, when Cartesian coordinates are used, *R* = 1 and *ξ* = 0, see FLOW3D manual FLOW3D (2007). *Ax* is the fractional area open to flow in the *x* direction, analogously for *Ay* and *Az*. Next, *VF* is the fractional volume open to flow, *R* and *ξ* are coefficients which depend on the coordinate system: (*x*, *y*, *z*) or (*r*, *θ*, *z*), *ρ* is the fluid density, *Rsor* is a mass source term. Finally, (*Gx*, *Gy*, *Gz*) are so called body accelerations FLOW3D (2007), (*fx*, *fy*, *fz*) are viscous 4 Will-be-set-by-IN-TECH

However, in the general case when there is no symmetry it is useful to apply the well known Cartesian coordinates *x*, *y*, and *z*: accordingly the equations of motion for the fluid velocity

In the specific case of the FLOW3D program the equations of motion for the fluid velocity components (*u*, *v*, *w*) with special additional terms included in the program are written:

> *Ayv*<sup>2</sup> *xVf*

*Ayuv xVf*

<sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂x*

<sup>=</sup> <sup>−</sup> *<sup>R</sup> ρ ∂p ∂y*

> *ρ ∂p ∂z*

*Rsor ρVf*

*Rsor ρVf*

> *Rsor ρVf*

*∂u ∂z* ) − *ξ*

*∂v ∂z* ) + *ξ*

+ *wAz*

Here, (*u*, *v*, *w*) are the velocity components in coordinate directions (*x*, *y*, *z*) respectively. For example, when Cartesian coordinates are used, *R* = 1 and *ξ* = 0, see FLOW3D manual FLOW3D (2007). *Ax* is the fractional area open to flow in the *x* direction, analogously for *Ay* and *Az*. Next, *VF* is the fractional volume open to flow, *R* and *ξ* are coefficients which depend on the coordinate system: (*x*, *y*, *z*) or (*r*, *θ*, *z*), *ρ* is the fluid density, *Rsor* is a mass source term. Finally, (*Gx*, *Gy*, *Gz*) are so called body accelerations FLOW3D (2007), (*fx*, *fy*, *fz*) are viscous

*∂w <sup>∂</sup><sup>z</sup>* ) = <sup>−</sup><sup>1</sup>

*∂vz*

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> 0, (7)

<sup>+</sup> *<sup>X</sup>* <sup>+</sup> *<sup>ν</sup>*∇2*u*, (8)

<sup>+</sup> *<sup>Y</sup>* <sup>+</sup> *<sup>ν</sup>*∇2*v*, (9)

<sup>+</sup> *<sup>Z</sup>* <sup>+</sup> *<sup>ν</sup>*∇2*w*. (10)

+ *Gx* + *fx* − *bx* −

+ *Gy* + *fy* − *by* −

+ *Gz* + *fz* − *bz* −

(*u* − *uw* − *δ* · *us*) (12)

(*v* − *vw* − *δ* · *vs*) (13)

(*w* − *ww* − *δ* · *ws*). (14)

*<sup>∂</sup><sup>z</sup>* <sup>=</sup> 0. (11)

(*rvr*) + <sup>1</sup> *r ∂v<sup>θ</sup> ∂θ* <sup>+</sup>

together with the continuity equation:

*∂u <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

*∂v <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

*∂w <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

In this case the continuity equation has the form:

components (*u*, *v*, *w*) are:

*∂u <sup>∂</sup><sup>t</sup>* <sup>+</sup>

*∂v <sup>∂</sup><sup>t</sup>* <sup>+</sup>

1 *VF* (*uAx ∂u ∂x*

1 *VF* (*uAx ∂v ∂x*

*∂w <sup>∂</sup><sup>t</sup>* <sup>+</sup>

1 *VF* (*uAx ∂w ∂x*

1 *r ∂ ∂r*

where *ν* = *η*/*ρ* is the kinematic viscosity Landau & Lifshitz (1959).

*∂u ∂x* + *v ∂u ∂y* + *w ∂u <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂x*

*∂v ∂x* + *v ∂v ∂y* + *w ∂v <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂y*

*∂w ∂x* + *v ∂w ∂y* + *w ∂w <sup>∂</sup><sup>z</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂z*

<sup>+</sup> *vAyR∂<sup>u</sup> ∂y*

<sup>+</sup> *vAyR∂<sup>v</sup> ∂y* *∂u ∂x* + *∂v ∂y* + *∂w*

+ *wAz*

+ *wAz*

<sup>+</sup> *vAyR∂<sup>w</sup> ∂y*

Fig. 1. The inflow time-dependent waveform used in these simulations. The function was taken from Fig. 3 of the work Papaharilaou et al. (2007).

accelerations, (*bx*, *by*, *bz*) are the flow losses in porous media or across porous baffle plates, and the final term accounts for the injection of mass at a source represented by a geometric component. Mathematical expressions for the viscous accelerations (*fx*, *fy*, *fz*) are presented in the Appendix.

The term *Uw* = (*uw*, *vw*, *ww*) in equations (12-14) is the velocity of the source component, which will generally be non-zero for a mass source of a General Moving Object (GMO) FLOW3D (2007). The term *Us* = (*us*, *vs*, *ws*) is the velocity of the fluid at the surface of the source relative to the source itself. It is computed in each control volume as

$$
\vec{\mathcal{U}}\_s = \frac{1}{\rho\_s} \frac{d(Q\vec{n})}{dA} \tag{15}
$$

where *dQ* is the mass flow rate, *ρ<sup>s</sup>* fluid source density, *dA* the area of the source surface in the cell and *n* the outward normal to the surface. The source is of the stagnation pressure type when in equations (12-14) *δ* = 0.0. Next, *δ* = 1.0 corresponds to the source of the static pressure type.

As we mentioned earlier, in all simulations we considered the blood flow as a pulsatile flow. The final result for the inflow waveform has been taken from Fig. 3 of the work Papaharilaou et al. (2007). The pulse was applied for 5.5 cycle times in our work. The waveform is shown in Fig. 1. These velocity values are used as time-dependent inflow initial boundary conditions. These numbers are included directly in the FLOW3D program.

Next, the general mass continuity equation, which is solved within the FLOW3D program has the following *general* form:

$$V\_f \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\rho u A\_x) + R \frac{\partial}{\partial y} (\rho v A\_y) + \frac{\partial}{\partial z} (\rho w A\_z) + \xi \frac{\rho u A\_x}{\mathbf{x}} = R\_{\text{dif}} + R\_{\text{sor}} \tag{16}$$

*Rdi f* is a turbulent diffusion term, and *Rsor* is a mass source. The turbulent diffusion term is

$$R\_{\rm dif} = \frac{\partial}{\partial \mathbf{x}} (v\_p A\_{\mathbf{x}} \frac{\partial \rho}{\partial \mathbf{x}}) + R \frac{\partial}{\partial y} (v\_p A\_y R \frac{\partial \rho}{\partial y}) + \frac{\partial}{\partial z} (v\_p A\_z \frac{\partial \rho}{\partial z}) + \xi \frac{\rho v\_p A\_{\mathbf{x}}}{\mathbf{x}},\tag{17}$$

components and etcetera. This procedure provides values for defining the flow parameters at discrete locations and allows specific boundary conditions to be set up. As the culminating step, one can start developing effective numerical approximations for the solution of the fluid

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 481

New pressure-velocity solutions have been implemented in FLOW-3D. We used the GMRES method. GMRES stands for the generalized minimum residual method. In addition to the GMRES solution, a new optional algorithm, the generalized conjugate gradient (GCG) algorithm, has also been implemented for solving viscous terms in the new GMRES routine. This new solver is a highly accurate and efficient method for a wide range of problems. It possesses good convergence, symmetry and speed properties; however, it does use more

Results of our simulations are presented below. One of the most important preliminary testing tasks is to check for numerical convergence. This test has been successfully accomplished in this work. A portion of the test calculation results are shown below in this paper. Next, in this work particular attention has been given to the calculations of the wall shear stress distribution (WSS). As we mentioned above WSS is the tangential drag force produced by moving blood, i.e. it is a mathematical function of the velocity gradient of blood near the

> *∂U*(*t*, *y*, *Rv*) *∂y*

Here *μ* is the dynamic viscosity, *t* is current time, *U*(*t*, *y*, *Rv*) is the flow velocity parallel to the wall, *y* is the distance to the wall of the vessel, and *Rv* is its radius. It was shown, that the magnitude of WSS is directly proportional to blood flow and blood viscosity and inversely proportional to the cube of the radius of the vessel, in other words a small change of the

First, we present results for a simple geometry vessel in the shape of a tube. However, the human blood is treated as real and a non-Newtonian liquid. The necessary data for viscosity

As we mentioned above, we take into account the real pulsatile flow, which is shown in Fig. 1. The data for Fig. 1 have also been obtained in clinical measurements, Papaharilaou et al. (2007). After such preliminary simulations we switched to a more complicated spatial configurations. This work begins with coronary bifurcation and the aortic arch. It is axiomatic that real people may have different size aortic arches with slightly different shapes. However,

The main goal of this work is to treat the above mentioned systems realistically, reveal the physics of the blood flow dynamics, and to obtain reliable results for pressure, dynamic viscosity, velocity profiles and strain rate distributions. Also, we tested the widely cited

First, we selected a simple vessel geometry, that is we considered the shape of a straight vessel to be a tube. In our simulations involving a straight cylinder type vessel we applied a cylindrical coordinate system: (*r*, *θ*, *Z*) with the axis *OZ* directed over the tube axis. Different quantities of cells have been used to discretize the empty space inside the tube. In the open space (inner part of the tube) the fluid dynamics equations have been solved using appropriate

*y*≈0

. (21)

*τ<sup>w</sup>* = *μ*

of the blood was found in previous laboratory and clinical measurements.

we carried out simulations for an average size and shape aortic arch.

Newtonian and non-Newtonian models of the human blood.

dynamics equations, i.e. NS equation.

memory than the SOR or SADI methods.

radius of a vessel will have a large effect on WSS.

**3. Numerical results**

endothelial surface:

**3.1 Straight vessel: cylinder**

where the coefficient *vp* = *Cpμ*/*ρ*, *μ* is dynamic viscosity and *Cp* is a constant. The *Rsor* term is a density source term that can be used to model mass injections through porous obstacle surfaces.

Compressible flow problems require the solution of the full density transport equation. In this work we treat blood as an incompressible fluid. For incompressible fluids *ρ* = *constant*, and the equation (15) becomes the following:

$$\frac{\partial}{\partial x}(uA\_x) + \frac{\partial}{\partial y}(vA\_y) + \frac{\partial}{\partial z}(wA\_z) + \xi \frac{uA\_x}{\varkappa} = \frac{R\_{\text{sor}}}{\rho}.\tag{18}$$

It is assumed, that at a stagnation pressure source fluid enters the domain at zero velocity. As a result, pressure should be considered at the source to move the fluid away from the source. For example, such sources are designed to model fluid emerging at the end of a rocket or the simple deflating process of a balloon. In general, stagnation pressure sources apply to cases when the momentum of the emerging fluid is created inside the source component, like in a rocket engine. At a static pressure source the fluid velocity is computed from the mass flow rate and the surface area of the source. In this case, no extra pressure is required to propel the fluid away from the source. An example of such a source is fluid emerging from a long straight pipe. Note that in this case the fluid momentum is created far from where the source is located.

Turbulence models can be taken into account in FLOW3D. It allows us to estimate the influence of turbulent fluctuations on mean flow quantities. This influence is usually expressed by additional diffusion terms in the equations for mean mass, momentum, and energy. The turbulence kinetic energy per unit mass, *q*, is the following:

$$\frac{\partial q}{\partial t} + \frac{1}{V\_F} \left( \mu A\_x \frac{\partial q}{\partial x} + vA\_y R \frac{\partial q}{\partial y} + wA\_z \frac{\partial q}{\partial z} \right) = P + G + Diff - D,\tag{19}$$

where *P* is shear production, *G* is buoyancy production, *Dif f* is diffusion, and *D* is a coefficient FLOW3D (2007).

When the turbulence option is used, the viscosity is a sum of the molecular and turbulent values. For non-Newtonian fluids the viscosity can be a function of the strain rate and/or temperature. A general expression based on the Carreau model is used in FLOW-3D for the strain rate dependent viscosity:

$$
\mu = \mu\_{\infty} + \frac{\mu\_0 E\_T - \mu\_{\infty}}{\lambda\_{00} + \left[\lambda\_0 + (\lambda\_1 E\_T)^2 e\_{ij} e\_{ij}\right]^{(1-n)/2}} + \frac{\lambda\_2}{\sqrt{(e\_{ij} e\_{ij})}},\tag{20}
$$

where *eij* = 1/2(*∂ui*/*∂xj* + *∂uj*/*∂xi*) is the fluid strain rate in Cartesian tensor notations, *μ*∞, *μ*0, *λ*0, *λ*1, *λ*<sup>2</sup> and *n* are constants. Also, *ET* = *exp*[*a*(*T*∗/(*T* − *b*) − *C*)], where *T*∗, *a*, *b*, and *C* are also parameters of the temperature dependence, and *T* is fluid temperature. This basic formula is used in our simulations for blood flow in vessels and in the aortic arch. For a variable dynamic viscosity *μ*, the viscous accelerations have a special form. That form is shown in the Appendix.

The equations of fluid dynamics should be solved together with specific boundary conditions. The numerical model starts with a computational mesh, or grid. It consists of a number of interconnected elements, or 3D-cells. These 3D-cells subdivide the physical space into small volumes with several nodes associated with each such volume. The nodes are used to store values of the unknown parameters, such as pressure, strain rate, temperature, velocity components and etcetera. This procedure provides values for defining the flow parameters at discrete locations and allows specific boundary conditions to be set up. As the culminating step, one can start developing effective numerical approximations for the solution of the fluid dynamics equations, i.e. NS equation.

New pressure-velocity solutions have been implemented in FLOW-3D. We used the GMRES method. GMRES stands for the generalized minimum residual method. In addition to the GMRES solution, a new optional algorithm, the generalized conjugate gradient (GCG) algorithm, has also been implemented for solving viscous terms in the new GMRES routine. This new solver is a highly accurate and efficient method for a wide range of problems. It possesses good convergence, symmetry and speed properties; however, it does use more memory than the SOR or SADI methods.

#### **3. Numerical results**

6 Will-be-set-by-IN-TECH

where the coefficient *vp* = *Cpμ*/*ρ*, *μ* is dynamic viscosity and *Cp* is a constant. The *Rsor* term is a density source term that can be used to model mass injections through porous obstacle

Compressible flow problems require the solution of the full density transport equation. In this work we treat blood as an incompressible fluid. For incompressible fluids *ρ* = *constant*, and

*∂z*

It is assumed, that at a stagnation pressure source fluid enters the domain at zero velocity. As a result, pressure should be considered at the source to move the fluid away from the source. For example, such sources are designed to model fluid emerging at the end of a rocket or the simple deflating process of a balloon. In general, stagnation pressure sources apply to cases when the momentum of the emerging fluid is created inside the source component, like in a rocket engine. At a static pressure source the fluid velocity is computed from the mass flow rate and the surface area of the source. In this case, no extra pressure is required to propel the fluid away from the source. An example of such a source is fluid emerging from a long straight pipe. Note that in this case the fluid momentum is created far from where the source

Turbulence models can be taken into account in FLOW3D. It allows us to estimate the influence of turbulent fluctuations on mean flow quantities. This influence is usually expressed by additional diffusion terms in the equations for mean mass, momentum, and

+ *wAz*

where *P* is shear production, *G* is buoyancy production, *Dif f* is diffusion, and *D* is a

When the turbulence option is used, the viscosity is a sum of the molecular and turbulent values. For non-Newtonian fluids the viscosity can be a function of the strain rate and/or temperature. A general expression based on the Carreau model is used in FLOW-3D for the

*<sup>λ</sup>*<sup>00</sup> + [*λ*<sup>0</sup> + (*λ*1*ET*)<sup>2</sup>*eijeij*](1−*n*)/2 <sup>+</sup>

where *eij* = 1/2(*∂ui*/*∂xj* + *∂uj*/*∂xi*) is the fluid strain rate in Cartesian tensor notations, *μ*∞, *μ*0, *λ*0, *λ*1, *λ*<sup>2</sup> and *n* are constants. Also, *ET* = *exp*[*a*(*T*∗/(*T* − *b*) − *C*)], where *T*∗, *a*, *b*, and *C* are also parameters of the temperature dependence, and *T* is fluid temperature. This basic formula is used in our simulations for blood flow in vessels and in the aortic arch. For a variable dynamic viscosity *μ*, the viscous accelerations have a special form. That form is

The equations of fluid dynamics should be solved together with specific boundary conditions. The numerical model starts with a computational mesh, or grid. It consists of a number of interconnected elements, or 3D-cells. These 3D-cells subdivide the physical space into small volumes with several nodes associated with each such volume. The nodes are used to store values of the unknown parameters, such as pressure, strain rate, temperature, velocity

*∂q ∂z* 

(*wAz*) + *ξ*

*uAx*

*<sup>x</sup>* <sup>=</sup> *Rsor*

*<sup>ρ</sup>* . (18)

= *P* + *G* + *Dif f* − *D*, (19)

, (20)

*λ*2 <sup>√</sup>(*eijeij*)

(*vAy*) + *<sup>∂</sup>*

surfaces.

is located.

*∂q <sup>∂</sup><sup>t</sup>* <sup>+</sup>

coefficient FLOW3D (2007).

strain rate dependent viscosity:

shown in the Appendix.

1 *VF*  *uAx ∂q ∂x*

the equation (15) becomes the following:

*∂ ∂x* (*uAx*) + *<sup>∂</sup>*

*∂y*

energy. The turbulence kinetic energy per unit mass, *q*, is the following:

*<sup>μ</sup>* <sup>=</sup> *<sup>μ</sup>*<sup>∞</sup> <sup>+</sup> *<sup>μ</sup>*0*ET* <sup>−</sup> *<sup>μ</sup>*<sup>∞</sup>

<sup>+</sup> *vAyR <sup>∂</sup><sup>q</sup> ∂y* Results of our simulations are presented below. One of the most important preliminary testing tasks is to check for numerical convergence. This test has been successfully accomplished in this work. A portion of the test calculation results are shown below in this paper. Next, in this work particular attention has been given to the calculations of the wall shear stress distribution (WSS). As we mentioned above WSS is the tangential drag force produced by moving blood, i.e. it is a mathematical function of the velocity gradient of blood near the endothelial surface:

$$\pi\_w = \mu \left[ \frac{\partial \mathcal{U}(t, y, R\_v)}{\partial y} \right]\_{y \approx 0}. \tag{21}$$

Here *μ* is the dynamic viscosity, *t* is current time, *U*(*t*, *y*, *Rv*) is the flow velocity parallel to the wall, *y* is the distance to the wall of the vessel, and *Rv* is its radius. It was shown, that the magnitude of WSS is directly proportional to blood flow and blood viscosity and inversely proportional to the cube of the radius of the vessel, in other words a small change of the radius of a vessel will have a large effect on WSS.

First, we present results for a simple geometry vessel in the shape of a tube. However, the human blood is treated as real and a non-Newtonian liquid. The necessary data for viscosity of the blood was found in previous laboratory and clinical measurements.

As we mentioned above, we take into account the real pulsatile flow, which is shown in Fig. 1. The data for Fig. 1 have also been obtained in clinical measurements, Papaharilaou et al. (2007). After such preliminary simulations we switched to a more complicated spatial configurations. This work begins with coronary bifurcation and the aortic arch. It is axiomatic that real people may have different size aortic arches with slightly different shapes. However, we carried out simulations for an average size and shape aortic arch.

The main goal of this work is to treat the above mentioned systems realistically, reveal the physics of the blood flow dynamics, and to obtain reliable results for pressure, dynamic viscosity, velocity profiles and strain rate distributions. Also, we tested the widely cited Newtonian and non-Newtonian models of the human blood.

#### **3.1 Straight vessel: cylinder**

First, we selected a simple vessel geometry, that is we considered the shape of a straight vessel to be a tube. In our simulations involving a straight cylinder type vessel we applied a cylindrical coordinate system: (*r*, *θ*, *Z*) with the axis *OZ* directed over the tube axis. Different quantities of cells have been used to discretize the empty space inside the tube. In the open space (inner part of the tube) the fluid dynamics equations have been solved using appropriate

0 0.2 0.4 0.6 0.8 1


**3.2 Hemodynamics in the coronary bifurcation**

be used.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Pulse time [sec]

Fig. 3. Time-dependent results for a specific geometrical point inside the cylinder: the middle point. Black dashed line: simulations without taking into account the turbulence; red bold line results with the turbulence. The non-Newtonian viscosity is taken into account.

significant, particularly in regard to dynamic viscosity and strain rate. This result means that in the case of pulsatile flows and non-Newtonian viscosity the turbulent term should be taken into account. In Fig. 4 we separately show the results for the pulsatile pressure distribution

Finally, it would also be very interesting to make a comparison between the results calculated using both a Newtonian and non-Newtonian viscosity. However, as in previous simulations, we will apply the pulsatile flow with the turbulence included, since it has proved to be important. The results are shown in Fig. 5. As one can see, we obtain significant differences between these two calculations. We specifically observed that for the pressure distribution,

Thus, we arrive at the important conclusion: within a time-dependent (pulsatile) flow of human blood it is necessary to take into account turbulence and non-Newtonian viscosity. The bold lines are results with the non-Newtonian viscosity (20) and the dashed lines are results with the Newtonian model when the viscosity *μ* has a constant value and is equal to 0.0345 P. As one can see the results are different for strain rate distributions and very different for pressure distributions. These results clearly indicate that in most cases when computer simulations are used in regard to human blood flow only the non-Newtonian model should

Below we show the result of a subsequent simulation involving a 90◦ bifurcated coronary artery in Figs. 6 and 7. The geometrical model of the bifurcation consisted of a 90◦ intersection of two cylinders. This model represents the bifurcation between the left anterior descending coronary artery and the circumflex coronary artery. In our opinion, in the case of pulsatile flow it is more interesting to present results in a time-dependent way. This method can provide a wider picture of highly non-stationary flow systems. In this paper, because of space limitations, we just included time-dependent results for pressure, dynamic viscosity,

and the turbulent energy, again using the middle point of the cylinder.

dynamic viscosity and turbulent energy, we obtained significant variation.

Velocity W [cm/sec]

Dynamic viscosity [poise]

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 483

Strain rate [1/sec]

Velocity V [cm/sec]

Fig. 2. Test of numerical convergence. Time-dependent dynamic viscosity, strain rate and velocity components V and W. Results for a vessel of simple geometry - cylinder type, for a specific spatial point inside the cylinder - the middle point. No turbulence effects are involved in these simulations with the realistic non-Newtonian viscosity of human blood. Black dashed line: calculations with 0.08 size for all cells FLOW3D (2007), red dot-dashed line with 0.07, green double dot - dashed line with 0.065, and blue bold line calculations with 0.062 size for all cells.

mathematical boundary conditions. The size of the tube is: *L* = 8 cm (in length) and *R* = 0.34 cm (length of inner radius). The thickness of the vessel wall is *s* = 0.03 cm. We have applied 5.5 cycles of blood pulse.

Let us now evaluate the expression (20). In these calculations we followed the work Cho & Kensey (1991), where the Carreau model of the human blood has also been used. To be consistent with Cho & Kensey (1991) we choose the following coefficients: *λ*<sup>2</sup> = *λ*<sup>00</sup> = 0, *a* = 0 and *ET* = 1, that is we don't take into account the temperature dependence of the viscosity. This investigation is to be addressed in our subsequent work. Next: *λ*<sup>0</sup> = 1, *λ*<sup>1</sup> = 3.313 sec, *μ*<sup>∞</sup> = 0.0345 P, *μ*<sup>0</sup> = 0.56 P, and *n* = 0.3568. The convergence was achieved when we used 52,800 cells, that is we used 100 points over *OZ*, 22 points over the radius of the inside space *R* = 0.34 cm, and 24 points over azimuthal angle Φ from 0 to 2*π*.

Time-dependent results for dynamic viscosity, strain rate and velocity components *V* and *W* are presented in Fig. 2. The turbulent effects are not taken into account. We decided to present only one precise geometrical point for comparison purposes: the middle point: *r* = *θ* = 0, and *Z* = 4.0 cm. The data for Fig. 2 were obtained with the non-Newtonian model of human blood. We refer the reader to the comments provided for the figure. We were able to closely replicate the values for all previous cell sizes FLOW3D (2007) and obtain almost identical values, for example for pressure, wall shear stress and other parameters, for 0.065 mm and 0.062 mm cell sizes FLOW3D (2007). This means, that the convergence has been achieved.

Next, it would be very interesting to compare the results calculated with and without the turbulent effect. To support this endeavor we used the realistic non-Newtonian model of blood viscosity, the pulsatile flow, and the size of computation cells at which convergence has been achieved, that is the 0.062 mm size for all computational cells FLOW3D (2007). The results are presented in Fig. 3. As we see from Fig. 3 the effect of the turbulence is

Fig. 3. Time-dependent results for a specific geometrical point inside the cylinder: the middle point. Black dashed line: simulations without taking into account the turbulence; red bold line results with the turbulence. The non-Newtonian viscosity is taken into account.

significant, particularly in regard to dynamic viscosity and strain rate. This result means that in the case of pulsatile flows and non-Newtonian viscosity the turbulent term should be taken into account. In Fig. 4 we separately show the results for the pulsatile pressure distribution and the turbulent energy, again using the middle point of the cylinder.

Finally, it would also be very interesting to make a comparison between the results calculated using both a Newtonian and non-Newtonian viscosity. However, as in previous simulations, we will apply the pulsatile flow with the turbulence included, since it has proved to be important. The results are shown in Fig. 5. As one can see, we obtain significant differences between these two calculations. We specifically observed that for the pressure distribution, dynamic viscosity and turbulent energy, we obtained significant variation.

Thus, we arrive at the important conclusion: within a time-dependent (pulsatile) flow of human blood it is necessary to take into account turbulence and non-Newtonian viscosity. The bold lines are results with the non-Newtonian viscosity (20) and the dashed lines are results with the Newtonian model when the viscosity *μ* has a constant value and is equal to 0.0345 P. As one can see the results are different for strain rate distributions and very different for pressure distributions. These results clearly indicate that in most cases when computer simulations are used in regard to human blood flow only the non-Newtonian model should be used.

#### **3.2 Hemodynamics in the coronary bifurcation**

8 Will-be-set-by-IN-TECH

Dynamic viscosity [poise]

Strain rate [1/sec]

Velocity, V [cm/sec]

Velocity, W [cm/sec]

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Pulse time [sec]

Fig. 2. Test of numerical convergence. Time-dependent dynamic viscosity, strain rate and velocity components V and W. Results for a vessel of simple geometry - cylinder type, for a specific spatial point inside the cylinder - the middle point. No turbulence effects are involved in these simulations with the realistic non-Newtonian viscosity of human blood. Black dashed line: calculations with 0.08 size for all cells FLOW3D (2007), red dot-dashed line with 0.07, green double dot - dashed line with 0.065, and blue bold line calculations with

mathematical boundary conditions. The size of the tube is: *L* = 8 cm (in length) and *R* = 0.34 cm (length of inner radius). The thickness of the vessel wall is *s* = 0.03 cm. We have applied

Let us now evaluate the expression (20). In these calculations we followed the work Cho & Kensey (1991), where the Carreau model of the human blood has also been used. To be consistent with Cho & Kensey (1991) we choose the following coefficients: *λ*<sup>2</sup> = *λ*<sup>00</sup> = 0, *a* = 0 and *ET* = 1, that is we don't take into account the temperature dependence of the viscosity. This investigation is to be addressed in our subsequent work. Next: *λ*<sup>0</sup> = 1, *λ*<sup>1</sup> = 3.313 sec, *μ*<sup>∞</sup> = 0.0345 P, *μ*<sup>0</sup> = 0.56 P, and *n* = 0.3568. The convergence was achieved when we used 52,800 cells, that is we used 100 points over *OZ*, 22 points over the radius of the inside space

Time-dependent results for dynamic viscosity, strain rate and velocity components *V* and *W* are presented in Fig. 2. The turbulent effects are not taken into account. We decided to present only one precise geometrical point for comparison purposes: the middle point: *r* = *θ* = 0, and *Z* = 4.0 cm. The data for Fig. 2 were obtained with the non-Newtonian model of human blood. We refer the reader to the comments provided for the figure. We were able to closely replicate the values for all previous cell sizes FLOW3D (2007) and obtain almost identical values, for example for pressure, wall shear stress and other parameters, for 0.065 mm and 0.062 mm cell sizes FLOW3D (2007). This means, that the convergence has been achieved. Next, it would be very interesting to compare the results calculated with and without the turbulent effect. To support this endeavor we used the realistic non-Newtonian model of blood viscosity, the pulsatile flow, and the size of computation cells at which convergence has been achieved, that is the 0.062 mm size for all computational cells FLOW3D (2007). The results are presented in Fig. 3. As we see from Fig. 3 the effect of the turbulence is

*R* = 0.34 cm, and 24 points over azimuthal angle Φ from 0 to 2*π*.

0 0.2 0.4 0.6

> -2 -1 0 1 2

0.062 size for all cells.

5.5 cycles of blood pulse.

Below we show the result of a subsequent simulation involving a 90◦ bifurcated coronary artery in Figs. 6 and 7. The geometrical model of the bifurcation consisted of a 90◦ intersection of two cylinders. This model represents the bifurcation between the left anterior descending coronary artery and the circumflex coronary artery. In our opinion, in the case of pulsatile flow it is more interesting to present results in a time-dependent way. This method can provide a wider picture of highly non-stationary flow systems. In this paper, because of space limitations, we just included time-dependent results for pressure, dynamic viscosity,

Presure [dynes/sq-cm]

0 0.2 0.4 0.6 0.8 1

0 3000 6000

> 0 0.2 0.4

approximation.

0 0.5 1 1.5 2 2.5 3 Pulse time [sec]

Time-dependent results for the middle point of the cylinder. Bold black line calculations with

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Pulse time [sec]

non-Newtonian viscosity, and the turbulence effect is included. Bold black line: results for the far right outflow side *z* = 0.0; red dashed line results for the farthest up outflow side

Once again we are using the Cartesian coordinate system. We also carried out a convergence test. To better represent the shape of the arch we applied five Cartesian sub-coordinate systems in our FLOW3D simulations. After the discretization the total number of all cubic cells reached about 900,000. It is important that once again we obtained full numerical convergence. Again, the geometry is shown in Figs. 8 and 9. In this work we computed pressure, velocity and strain rate distributions in the arch, while the human blood is treated

Fig. 5. Results for pressure, dynamic viscosity, turbulent energy and velocity W.

Pressure [dynes/sq-cm]

Dynamic viscosity [poise]

Turbulent energy [ergs/gm]

Fig. 6. Time-dependent results for a vessel with bifurcation. Pulsatile blood flow,

as a non-Newtonian liquid and while the realistic pulsatile blood flow is used.

Strain rate [1/sec]

*y* = 0.0. Convergence test results for velocity components: U, V, W.

non-Newtonian viscosity of the human blood; red dashed line with its Newtonian

Dynamic viscosity [poise]

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 485

Velocity W [cm/sec]

Turbulent energy [ergs/gm]

Fig. 4. Time-dependent results for pressure and the turbulent energy in the middle point of the cylinder. The non-Newtonian viscosity.

turbulent energy, and strain rate. However, we understand, that results which depend on spatial coordinates (*r*, *θ*, *Z*) for a few fixed moments of time are also highly useful.

In the case of the bifurcation shown in Fig. 6, we report the results for only two spatial points, which are the two outflow sides: the far right side and the farthest upper side of the bifurcation. The length of the lower horizontal vessel is 4 cm and its diameter is 0.54 cm. The length of the upper vertical vessel is 1.2 cm and its diameter is 0.4 cm. These sizes are consistent with average size human vessels.

Further, Fig. 6 represents our time-dependent results for the two outflow sides mentioned above. These results are for pressure, dynamic viscosity, turbulent energy and strain rate. The bold black lines are the results for the right outflow side, and the red dashed lines are the results for the farthest upper side (see comments to Fig. 6). In conclusion, the main goal of these calculations is to adopt them to investigate a case in which a stent is implanted in the bifurcation area Frank et al. (2002).

In Fig. 7 blood flows in from the left to the right with the imposed initial velocity profile taken from Fig. 1. The pressure, strain rate and turbulent energy distributions are shown for only one specific time moment *t*=4.329 s. The velocity vectors are also shown on these plots.

#### **3.3 Blood flow in aortic arch**

The geometry of the blood simulations inside the human aortic arch is shown in Figs. 8 and 9. On the top of the aortic arch three arteries are included. These arteries deliver the blood to the carotid artery and then to the brain. This configuration only models and approximately represents the real aortic arch. One of the goals of our simulations is to reveal the physics of the blood flow dynamics in this important portion of the human cardiovascular system.

The aortic arch is represented as a curved tube. The outer radius of the tube is 2.6 cm. A straight vessel (tube) is also merged to the arch. The length of the straight tube is about 4 cm. Again, the thickness of the wall is 0.03 cm, and the inner radius of the tube is *r* = 0.34 cm. The thickness is not important in these simulations, but it will be useful when, in future works, we will need to introduce elasticity of the walls of the tubes. The FLOW3D program allows to carry out fluid dynamic simulations with elastic (not only hard body) walls.

10 Will-be-set-by-IN-TECH

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 Pulse time [sec]

Fig. 4. Time-dependent results for pressure and the turbulent energy in the middle point of

turbulent energy, and strain rate. However, we understand, that results which depend on

In the case of the bifurcation shown in Fig. 6, we report the results for only two spatial points, which are the two outflow sides: the far right side and the farthest upper side of the bifurcation. The length of the lower horizontal vessel is 4 cm and its diameter is 0.54 cm. The length of the upper vertical vessel is 1.2 cm and its diameter is 0.4 cm. These sizes are

Further, Fig. 6 represents our time-dependent results for the two outflow sides mentioned above. These results are for pressure, dynamic viscosity, turbulent energy and strain rate. The bold black lines are the results for the right outflow side, and the red dashed lines are the results for the farthest upper side (see comments to Fig. 6). In conclusion, the main goal of these calculations is to adopt them to investigate a case in which a stent is implanted in the

In Fig. 7 blood flows in from the left to the right with the imposed initial velocity profile taken from Fig. 1. The pressure, strain rate and turbulent energy distributions are shown for only one specific time moment *t*=4.329 s. The velocity vectors are also shown on these plots.

The geometry of the blood simulations inside the human aortic arch is shown in Figs. 8 and 9. On the top of the aortic arch three arteries are included. These arteries deliver the blood to the carotid artery and then to the brain. This configuration only models and approximately represents the real aortic arch. One of the goals of our simulations is to reveal the physics of the blood flow dynamics in this important portion of the human cardiovascular system. The aortic arch is represented as a curved tube. The outer radius of the tube is 2.6 cm. A straight vessel (tube) is also merged to the arch. The length of the straight tube is about 4 cm. Again, the thickness of the wall is 0.03 cm, and the inner radius of the tube is *r* = 0.34 cm. The thickness is not important in these simulations, but it will be useful when, in future works, we will need to introduce elasticity of the walls of the tubes. The FLOW3D program allows to

carry out fluid dynamic simulations with elastic (not only hard body) walls.

spatial coordinates (*r*, *θ*, *Z*) for a few fixed moments of time are also highly useful.

the cylinder. The non-Newtonian viscosity.

consistent with average size human vessels.

bifurcation area Frank et al. (2002).

**3.3 Blood flow in aortic arch**

Turbulent energy [ergs/gm]

Pressure [dynes/cm2

]

Fig. 5. Results for pressure, dynamic viscosity, turbulent energy and velocity W. Time-dependent results for the middle point of the cylinder. Bold black line calculations with non-Newtonian viscosity of the human blood; red dashed line with its Newtonian approximation.

Fig. 6. Time-dependent results for a vessel with bifurcation. Pulsatile blood flow, non-Newtonian viscosity, and the turbulence effect is included. Bold black line: results for the far right outflow side *z* = 0.0; red dashed line results for the farthest up outflow side *y* = 0.0. Convergence test results for velocity components: U, V, W.

Once again we are using the Cartesian coordinate system. We also carried out a convergence test. To better represent the shape of the arch we applied five Cartesian sub-coordinate systems in our FLOW3D simulations. After the discretization the total number of all cubic cells reached about 900,000. It is important that once again we obtained full numerical convergence. Again, the geometry is shown in Figs. 8 and 9. In this work we computed pressure, velocity and strain rate distributions in the arch, while the human blood is treated as a non-Newtonian liquid and while the realistic pulsatile blood flow is used.

0

0

167

335

503

671

 -4.0 -2.5 -1.0 0.5 2.0 3.5 x \*cm\*

strain rate (1/s)

 -4.0 -2.5 -1.0 0.5 2.0 3.5 x \*cm\*

Fig. 8. Blood flow in the aortic arch. These two plots represent the full 2D-picture of the geometry used in these simulations. Shaded results for the strain rate are also shown, the bars on the right show the values. Results are for two specific moments of the time *t*<sup>40</sup> = 4.329 sec and *t*<sup>41</sup> = 4.440 sec. The values of the strain rate distribution range from 0.0 1/sec to 357.0 1/sec (upper plot) and from 0.0 to 671 1/sec (lower plot). The maximum values of the strain rate are localized in the region inside the arch. Blood flows from right to left in both pictures.

strain rate (1/s)

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 487


Title

2.6

y \* c m \*


Title

m-b linked

m-b linked

FLOW-3D t=4.329 z=2.500E-02 ix=3 to 20 jy=2 to 43

FLOW-3D t=4.440 z=2.500E-02 ix=3 to 20 jy=2 to 43

19:56:38 06/05/2008 jbyh hydr3d: version 9.2 win32-ifl 2006

19:56:38 06/05/2008 jbyh hydr3d: version 9.2 win32-ifl 2006

2.6

y \* c m \*

89

178

267

357

Fig. 7. The figures are 2D-plots showing the blood flow in the bifurcated vessels for only one precise moment of the discretized time *ti* = 4.329 sec, the corresponding index is *i* = 40. Upper plot represents the result for the pressure distribution in the bifurcation, and the pressure ranges from 2068 dynes/sq-cm to 6758 dynes/sq-cm. The middle plot represents the results for the strain rate distribution and the lower plot shows results for the turbulent energy in the bifurcation. The range of the values is also shown.

12 Will-be-set-by-IN-TECH

pressure (dynes/sq-cm)

 -4.0 -3.2 -2.4 -1.6 -0.8 0.0 x \*cm\*

strain rate (1/s)

 -4.0 -3.2 -2.4 -1.6 -0.8 0.0 x \*cm\*

turbulent energy (ergs/gm)

 -4.0 -3.2 -2.4 -1.6 -0.8 0.0 x \*cm\*

Fig. 7. The figures are 2D-plots showing the blood flow in the bifurcated vessels for only one precise moment of the discretized time *ti* = 4.329 sec, the corresponding index is *i* = 40. Upper plot represents the result for the pressure distribution in the bifurcation, and the pressure ranges from 2068 dynes/sq-cm to 6758 dynes/sq-cm. The middle plot represents the results for the strain rate distribution and the lower plot shows results for the turbulent


Title


Title



Title

m-b linked

m-b linked


m-b linked

FLOW-3D t=4.329 z=2.700E-02 ix=2 to 9 jy=3 to 25

FLOW-3D t=4.329 z=2.700E-02 ix=2 to 9 jy=3 to 25

FLOW-3D t=4.329 z=2.700E-02 ix=2 to 9 jy=3 to 25

energy in the bifurcation. The range of the values is also shown.

14:46:35 06/11/2008 nomj hydr3d: version 9.2 win32-ifl 2006

14:46:35 06/11/2008 nomj hydr3d: version 9.2 win32-ifl 2006

14:46:35 06/11/2008 nomj hydr3d: version 9.2 win32-ifl 2006


2068

1

0

754

1508

2262

3016

(max=4.67E+01)

96

190

285

380

3240

4413

5585

6758

(max=4.67E+01)

(max=4.67E+01)

Fig. 8. Blood flow in the aortic arch. These two plots represent the full 2D-picture of the geometry used in these simulations. Shaded results for the strain rate are also shown, the bars on the right show the values. Results are for two specific moments of the time *t*<sup>40</sup> = 4.329 sec and *t*<sup>41</sup> = 4.440 sec. The values of the strain rate distribution range from 0.0 1/sec to 357.0 1/sec (upper plot) and from 0.0 to 671 1/sec (lower plot). The maximum values of the strain rate are localized in the region inside the arch. Blood flows from right to left in both pictures.

In Fig. 8 we present the results of strain rate distributions inside the arch for two specific time moments. At the most left point, which is the inlet, we specify the pulsatile velocity source as the initial condition, that is the data from Fig. 1 are used. From the general theory of fluid mechanics Landau & Lifshitz (1959) it is possible to determine together with the blood density, viscosity, and spatial geometries, the dynamics of the blood according to the Navier-Stokes equation and its boundary conditions. Small vectors indicate the blood velocity. As can be seen from Fig. 8 blood flows from left to right in direction. However, because of pulsatility

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 489

The values of the strain rate are also shown. These values are strongly oscillating. From the plots one can conclude that in the region of the arch the strain rate values become much larger than in the region of the straight vessel. This result represents clear evidence that in this part of the human vascular system atherosclerotic plaques should localize less than in the straight vessels. However, the higher wall shear stress values in the aortic arch could be the reason for sudden mechanical disruption of the arterial wall in this part of the human vascular system. These results are consistent with laboratory and clinical observations. In Fig. 9 we depict the

In this work we applied computational fluid dynamics techniques to support pulsatile human blood flow simulations through different shape/size vessels and the aortic arch. The realistic blood pulse has been adopted and applied from the work Papaharilaou et al. (2007). The geometrical size of the vessels and the aortic arch have been selected to match the average real values. Human blood was treated in two different ways: (a) as a Newtonian liquid when the viscosity of the blood has a constant value, and (b) as a non-Newtonian liquid with the viscosity value represented by the equation (20). The numerical coefficients in (20) have been

It is always difficult to obtain a steady-state cycle profile and stable computational results at the very beginning of time-dependent simulations. However, after a short stabilization period a steady-state cycle profile can be obtained. In our simulations we used up to 5.5 pulse cycles to reach complete steady state profiles. We obtained valid results for pressure, wall shear stress distribution and other physical parameters, such as the three velocity components of

Our simulations showed that the FLOW3D program is capable of providing stable numerical results for all geometries included in this work. The time-dependent mathematical convergence test has been successfully carried out. Particular attention has been paid to this aspect of the calculations. It is a well known fact that fluid dynamics equations can have unstable solutions, Landau & Lifshitz (1959). Therefore, numerical convergence has been

The result of computer simulations of blood flow in vessels for three different geometries have been presented. For pressure, strain rate and velocity component distributions we found significant disagreements between our results obtained with the realistic non-Newtonian treatment of human blood and the widely used method in literature: a simple Newtonian

Our results are in good agreement with the conclusions of the works, Chen & Lu (2004; 2006), where the authors also obtained significant differences between their results calculated with and without the non-Newtonian effect of blood viscosity. However, the recent work, Boyd & Buick (2007) should be mentioned, in which the authors performed 2-dimensional simulations of human blood flow through the carotid artery with and without the non-Newtonian effect of

blood flows in the opposite direction too.

pressure distribution in the arch.

taken from work, Cho & Kensey (1991).

tested and confirmed in this work.

approximation.

blood flow. All of these were shown in Figs. 2-6.

**4. Conclusion**

Fig. 9. Blood flow in the aortic arch. These two plots represent in more detail the region of the arch together with shaded results for the pressure distribution. The bars on the right show the values. These results are for two specific moments of the time *t*<sup>40</sup> = 4.329 sec and *t*<sup>41</sup> = 4.440 sec, where the pressure ranges from 957 dynes/sp-cm to 5889 dynes/sq-cm (upper plot), and from -616 dynes/sq-cm to 7070 dynes/sq-cm (lower plot).

In Fig. 8 we present the results of strain rate distributions inside the arch for two specific time moments. At the most left point, which is the inlet, we specify the pulsatile velocity source as the initial condition, that is the data from Fig. 1 are used. From the general theory of fluid mechanics Landau & Lifshitz (1959) it is possible to determine together with the blood density, viscosity, and spatial geometries, the dynamics of the blood according to the Navier-Stokes equation and its boundary conditions. Small vectors indicate the blood velocity. As can be seen from Fig. 8 blood flows from left to right in direction. However, because of pulsatility blood flows in the opposite direction too.

The values of the strain rate are also shown. These values are strongly oscillating. From the plots one can conclude that in the region of the arch the strain rate values become much larger than in the region of the straight vessel. This result represents clear evidence that in this part of the human vascular system atherosclerotic plaques should localize less than in the straight vessels. However, the higher wall shear stress values in the aortic arch could be the reason for sudden mechanical disruption of the arterial wall in this part of the human vascular system. These results are consistent with laboratory and clinical observations. In Fig. 9 we depict the pressure distribution in the arch.

#### **4. Conclusion**

14 Will-be-set-by-IN-TECH

pressure (dynes/sq-cm)

 0.0 0.52 1.04 1.56 2.08 2.60 x \*cm\*

pressure (dynes/sq-cm)

 0.0 0.52 1.04 1.56 2.08 2.60 x \*cm\*

Fig. 9. Blood flow in the aortic arch. These two plots represent in more detail the region of the arch together with shaded results for the pressure distribution. The bars on the right show the values. These results are for two specific moments of the time *t*<sup>40</sup> = 4.329 sec and *t*<sup>41</sup> = 4.440 sec, where the pressure ranges from 957 dynes/sp-cm to 5889 dynes/sq-cm (upper


Title

1.40

y \* c m \*


Title

mb 2 linked

FLOW-3D t=4.440 z=2.500E-02 ix=2 to 53 jy=3 to 58

19:56:38 06/05/2008 jbyh hydr3d: version 9.2 win32-ifl 2006

plot), and from -616 dynes/sq-cm to 7070 dynes/sq-cm (lower plot).

mb 2 linked

FLOW-3D t=4.329 z=2.500E-02 ix=2 to 53 jy=3 to 58

19:56:38 06/05/2008 jbyh hydr3d: version 9.2 win32-ifl 2006

1.40

y \* c m \*

957


664

1945

3226

4507

5789

7070

(max=7.94E+01)

1779

2601

3423

4245

5067

5889

(max=4.27E+01)

In this work we applied computational fluid dynamics techniques to support pulsatile human blood flow simulations through different shape/size vessels and the aortic arch. The realistic blood pulse has been adopted and applied from the work Papaharilaou et al. (2007). The geometrical size of the vessels and the aortic arch have been selected to match the average real values. Human blood was treated in two different ways: (a) as a Newtonian liquid when the viscosity of the blood has a constant value, and (b) as a non-Newtonian liquid with the viscosity value represented by the equation (20). The numerical coefficients in (20) have been taken from work, Cho & Kensey (1991).

It is always difficult to obtain a steady-state cycle profile and stable computational results at the very beginning of time-dependent simulations. However, after a short stabilization period a steady-state cycle profile can be obtained. In our simulations we used up to 5.5 pulse cycles to reach complete steady state profiles. We obtained valid results for pressure, wall shear stress distribution and other physical parameters, such as the three velocity components of blood flow. All of these were shown in Figs. 2-6.

Our simulations showed that the FLOW3D program is capable of providing stable numerical results for all geometries included in this work. The time-dependent mathematical convergence test has been successfully carried out. Particular attention has been paid to this aspect of the calculations. It is a well known fact that fluid dynamics equations can have unstable solutions, Landau & Lifshitz (1959). Therefore, numerical convergence has been tested and confirmed in this work.

The result of computer simulations of blood flow in vessels for three different geometries have been presented. For pressure, strain rate and velocity component distributions we found significant disagreements between our results obtained with the realistic non-Newtonian treatment of human blood and the widely used method in literature: a simple Newtonian approximation.

Our results are in good agreement with the conclusions of the works, Chen & Lu (2004; 2006), where the authors also obtained significant differences between their results calculated with and without the non-Newtonian effect of blood viscosity. However, the recent work, Boyd & Buick (2007) should be mentioned, in which the authors performed 2-dimensional simulations of human blood flow through the carotid artery with and without the non-Newtonian effect of

*τxy* = −*μ*

*τxz* = −*μ*

*τyz* = −*μ*

fractional flow areas (*Ax*, *Ay*, *Az*) which vanish at the walls FLOW3D (2007).

In the above equations (22)-(24) the terms *w<sup>s</sup>*

Engineering 34 (8), pp. 1259-1271.

Springer-Verlag New-York, LLC.

Biomedical Materials 2 (1), art. no. S05, S28-S37.

with Elastic Boundaries, Phys. Rev. E 57(1), R25-R28.

shear stress magnitude.

1899-1911.

818-832.

pp. 357-382.

**7. References**

 *∂v <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>R</sup> <sup>∂</sup><sup>u</sup>*

Numerical Modeling and Simulations of Pulsatile Human Blood Flow in Different 3D-Geometries 491

 *∂u <sup>∂</sup><sup>z</sup>* <sup>+</sup> *∂w ∂x* 

 *∂v <sup>∂</sup><sup>z</sup>* <sup>+</sup> *<sup>R</sup> <sup>∂</sup><sup>w</sup> ∂y* 

*<sup>x</sup>*, *w<sup>s</sup>*

are equal to zero, there is no wall shear stress. This is because the remaining terms contain the

The wall stresses are modeled by assuming a zero tangential velocity on the portion of any area closed to flow. Mesh and moving obstacle boundaries are an exception because they can be assigned non-zero tangential velocities. In this case the allowed boundary motion corresponds to a rigid body translation of the boundary parallel to its surface. For turbulent flows, a law-of-the-wall velocity profile is assumed near the wall, which modifies the wall

Agarwal, R., Katiyar, V.K. & Pradhan, P. (2008). A mathematical modeling of pulsatile flow in carotid artery bifurcation, Intern. J. of Engineering Science 46 (11), pp. 1147-1156. Benard, N., Perrault, R. & Coisne, D. (2006). Computational Approach to Estimating the

Banerjee, R.K., Devarakonda, S.B., Rajamohan, D. & Back, L.H. (2007). Developed Pulsatile

Boyd, J. & Buick, J.M. (2007). Comparison of Newtonian and non-Newtonian Flows in

Carter, Y.M, Karmy-Jones, R.C., Oxorn, D.C. & Aldea, G.S. (2001).Traumatic Disruption of the

Chen, J. & Lu, X.-Y. (2005). Numerical investigation of the non-Newtonian blood flow in a

Chen, J. & Lu, X.-Y. (2006). Numerical investigation of the non-Newtonian pulsatile blood flow

Cho, Y.I. & Kensey, K.R. (1991). Effects of the non-Newtonian viscosity of blood on flows in a diseased arterial vessel. Part 1: Steady Flows, Biorheology 28(3-4), pp. 241-262. Dhein, S., Delmar, M., & Mohr, F.W. (2005). Practical Methods in Cardiovascular Research,

Duraiswamy, N., Schoephoerster, R.T., Moreno, M.R. & Moore Jr., J.E. (2007). Stented Artery

Faik, I., Mongrain, R., Leask, R.L., Rodes-Cabau, J., Larose, E. & Bertrand, O. (2007).

Fang, H., Lin,Z. & Wang,Z. (1998). Lattice Boltzmann Simulation of Viscous Fluid Systems

FLOW-3D Users Manual. (2007). Version 9.2, Flow Science, Santa Fe, New Mexico, USA.

Flow in a Deployed Coronary Stent, Biorheology 44 (2), pp. 91-102.

Physics in Medicine and Biology 52(20), pp. 6215-6228.

Aortic Arch, European J. Cardiothoracic Surgery, 20, pp. 1231.

Effects of Blood Properties on Changes in Intra-Stent Flow, Annals of Biomed.

a Two-Dimensional Carotid Artery Model Using the Lattice Boltzmann Method,

bifurcation model with a non-planar branch Journal of Biomechanics, 37 (12), pp.

in a bifurcation model with a non-planar branch, Journal of Biomechanics 39 (5), pp.

Flow Patterns and Their Effects on the Artery Wall, Ann. Rev. of Fluid Mechanics 39,

Time-Dependent 3D Simulations of the Hemodynamics in a Stented Coronary Artery,

*<sup>y</sup>* and *w<sup>s</sup>*

*<sup>∂</sup><sup>y</sup>* <sup>−</sup> *<sup>ξ</sup><sup>v</sup> x* 

(28)

(29)

. (30)

*<sup>z</sup>* are wall shear stresses. If these terms

the viscosity. They did not find any substantial differences in their results. Finally, we would like to mention the paper Agarwal et. al. (2008), where the authors also performed simulations for the carotid artery, but only the non-Newtonian viscosity was used.

Next, the influence of a possible turbulent effect has also been investigated in this work. It was found that the effect is important. We believe, that the physical reason of this phenomena lies in the strong pulsatility of the flow and in the non-Newtonian viscosity of the blood. The contribution of the turbulence is most significant in the area of bifurcated vessels.

Finally, a significant increase of the strain rate and, the wall shear stress distribution, is found in the region of the aortic arch. This computational result provides additional evidence to support recent clinical and laboratory observations that this part of the human cardiovascular system is under higher risk of disruption Carter et al. (2001); Pochettino & Bavaria (2006). In future work it would be interesting to include the elasticity of the walls of the aortic arch Fang et al. (1998) and other vessels.

In conclusion, we would like to specifically point out, that the developments in this work can be directly applied to even more interesting and very important situations such as when a stent is implanted inside a vessel Frank et al. (2002). In this case, for example, it would be very useful to determine blood flow disturbance, the pressure distribution, strain rate and values of other physical parameters. The results of this work should allow us to determine the optimal size and shape of effective stents. As we mentioned in the Introduction some research groups are carrying out laboratory and computer simulations of blood flow through vessels with implanted stents Frank et al. (2002). It is very difficult to underestimate the value of these works.

#### **5. Acknowledgments**

This work was partially supported by Office of Sponsored Programs (OSP), by Internal Grant Program of St. Cloud State University, St. Cloud, MN-56301-4498, USA, and by a private Minnesota based company: *RIE Coatings, Eden Valley*, MN-55329-1646, USA (www.riecoatings.com).

#### **6. Appendix**

For a variable dynamic viscosity *μ*, the viscous accelerations are

$$\rho V\_{\rm F} f\_{\rm x} = w\_{\rm x}^{\rm s} - \left[ \frac{\partial}{\partial \mathbf{x}} (A\_{\rm x} \tau\_{\rm xx}) + R \frac{\partial}{\partial y} (A\_{\rm y} \tau\_{\rm xy}) + \frac{\partial}{\partial z} (A\_{\rm z} \tau\_{\rm xz}) + \frac{\tilde{\xi}}{\mathbf{x}} (A\_{\rm x} \tau\_{\rm xx} - A\_{\rm y} \tau\_{\rm yy}) \right] \tag{22}$$

$$\rho V\_{\rm F} f\_{\rm y} = w\_{\rm y}^{s} - \left[ \frac{\partial}{\partial \mathbf{x}} (A\_{\rm x} \pi\_{\rm xy}) + R \frac{\partial}{\partial y} (A\_{\rm y} \pi\_{\rm yy}) + \frac{\partial}{\partial \mathbf{z}} (A\_{\rm z} \pi\_{\rm yz}) + \frac{\xi}{\mathbf{x}} (A\_{\rm x} + A\_{\rm y} \pi\_{\rm xy}) \right] \tag{23}$$

$$
\rho V\_{\rm F} f\_{\rm z} = w\_{\rm z}^s - \left[ \frac{\partial}{\partial x} (A\_{\rm x} \pi\_{\rm x}) + R \frac{\partial}{\partial y} (A\_{\rm y} \pi\_{\rm y}) + \frac{\partial}{\partial z} (A\_{\rm z} \pi\_{\rm z}) + \frac{\tilde{\xi}}{\mathfrak{x}} (A\_{\rm x} \pi\_{\rm x}) \right] \tag{24}
$$

where

$$\tau\_{\text{xx}} = -2\mu \left( \frac{\partial u}{\partial x} - \frac{1}{3} \left( \frac{\partial u}{\partial x} + R \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} + \frac{\xi u}{x} \right) \right) \tag{25}$$

$$\tau\_{yy} = -2\mu \left[ R \frac{\partial v}{\partial x} + \xi \frac{u}{x} - \frac{1}{3} \left( \frac{\partial u}{\partial x} + R \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} + \frac{\xi u}{x} \right) \right] \tag{26}$$

$$\tau\_{zz} = -2\mu \left( \frac{\partial w}{\partial z} - \frac{1}{3} \left( \frac{\partial u}{\partial x} + R \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} + \frac{\xi u}{x} \right) \right) \tag{27}$$

$$\tau\_{xy} = -\mu \left( \frac{\partial v}{\partial x} + R \frac{\partial u}{\partial y} - \frac{\xi v}{x} \right) \tag{28}$$

$$
\pi\_{\rm xz} = -\mu \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \tag{29}
$$

$$
\pi\_{yz} = -\mu \left( \frac{\partial v}{\partial z} + R \frac{\partial w}{\partial y} \right). \tag{30}
$$

In the above equations (22)-(24) the terms *w<sup>s</sup> <sup>x</sup>*, *w<sup>s</sup> <sup>y</sup>* and *w<sup>s</sup> <sup>z</sup>* are wall shear stresses. If these terms are equal to zero, there is no wall shear stress. This is because the remaining terms contain the fractional flow areas (*Ax*, *Ay*, *Az*) which vanish at the walls FLOW3D (2007).

The wall stresses are modeled by assuming a zero tangential velocity on the portion of any area closed to flow. Mesh and moving obstacle boundaries are an exception because they can be assigned non-zero tangential velocities. In this case the allowed boundary motion corresponds to a rigid body translation of the boundary parallel to its surface. For turbulent flows, a law-of-the-wall velocity profile is assumed near the wall, which modifies the wall shear stress magnitude.

#### **7. References**

16 Will-be-set-by-IN-TECH

the viscosity. They did not find any substantial differences in their results. Finally, we would like to mention the paper Agarwal et. al. (2008), where the authors also performed simulations

Next, the influence of a possible turbulent effect has also been investigated in this work. It was found that the effect is important. We believe, that the physical reason of this phenomena lies in the strong pulsatility of the flow and in the non-Newtonian viscosity of the blood. The

Finally, a significant increase of the strain rate and, the wall shear stress distribution, is found in the region of the aortic arch. This computational result provides additional evidence to support recent clinical and laboratory observations that this part of the human cardiovascular system is under higher risk of disruption Carter et al. (2001); Pochettino & Bavaria (2006). In future work it would be interesting to include the elasticity of the walls of the aortic arch Fang

In conclusion, we would like to specifically point out, that the developments in this work can be directly applied to even more interesting and very important situations such as when a stent is implanted inside a vessel Frank et al. (2002). In this case, for example, it would be very useful to determine blood flow disturbance, the pressure distribution, strain rate and values of other physical parameters. The results of this work should allow us to determine the optimal size and shape of effective stents. As we mentioned in the Introduction some research groups are carrying out laboratory and computer simulations of blood flow through vessels with implanted stents Frank et al. (2002). It is very difficult to underestimate the value of these

This work was partially supported by Office of Sponsored Programs (OSP), by Internal Grant Program of St. Cloud State University, St. Cloud, MN-56301-4498, USA, and by a private Minnesota based company: *RIE Coatings, Eden Valley*, MN-55329-1646, USA

(*Ayτxy*) + *<sup>∂</sup>*

(*Ayτyy*) + *<sup>∂</sup>*

(*Ayτyz*) + *<sup>∂</sup>*

*∂z*

*∂z*

*∂z*

(*Azτxz*) + *<sup>ξ</sup>*

(*Azτyz*) + *<sup>ξ</sup>*

(*Azτzz*) + *<sup>ξ</sup>*

*x*

*x*

*x*

(*Axτxx* − *Ayτyy*)] (22)

(*Ax* + *Ayτxy*)] (23)

(*Axτxz*)], (24)

(25)

(26)

(27)

For a variable dynamic viscosity *μ*, the viscous accelerations are

(*Axτxx*) + *<sup>R</sup> <sup>∂</sup>*

(*Axτxy*) + *<sup>R</sup> <sup>∂</sup>*

(*Axτxz*) + *<sup>R</sup> <sup>∂</sup>*

 *∂u <sup>∂</sup><sup>x</sup>* <sup>−</sup> <sup>1</sup> 3 *∂u <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>R</sup> <sup>∂</sup><sup>v</sup> ∂y* + *∂w <sup>∂</sup><sup>z</sup>* <sup>+</sup> *<sup>ξ</sup><sup>u</sup> x* 

 *∂w <sup>∂</sup><sup>z</sup>* <sup>−</sup> <sup>1</sup> 3 *∂u <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>R</sup> <sup>∂</sup><sup>v</sup> ∂y* + *∂w <sup>∂</sup><sup>z</sup>* <sup>+</sup> *<sup>ξ</sup><sup>u</sup> x* 

 *<sup>R</sup> <sup>∂</sup><sup>v</sup> <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>ξ</sup> u <sup>x</sup>* <sup>−</sup> <sup>1</sup> 3 *∂u <sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>R</sup> <sup>∂</sup><sup>v</sup> ∂y* + *∂w <sup>∂</sup><sup>z</sup>* <sup>+</sup> *<sup>ξ</sup><sup>u</sup> x* 

*∂y*

*∂y*

*∂y*

contribution of the turbulence is most significant in the area of bifurcated vessels.

for the carotid artery, but only the non-Newtonian viscosity was used.

et al. (1998) and other vessels.

**5. Acknowledgments**

(www.riecoatings.com).

*ρVF fx* = *w<sup>s</sup>*

*ρVF fy* = *w<sup>s</sup>*

*ρVF fz* = *w<sup>s</sup>*

*<sup>x</sup>* <sup>−</sup> [ *<sup>∂</sup> ∂x*

*<sup>y</sup>* <sup>−</sup> [ *<sup>∂</sup> ∂x*

*<sup>z</sup>* <sup>−</sup> [ *<sup>∂</sup> ∂x*

*τxx* = −2*μ*

*τyy* = −2*μ*

*τzz* = −2*μ*

**6. Appendix**

where

works.


**22** 

*Brazil* 

**Biomechanical Factors** 

Kleiber Bessa1, Daniel Legendre2 and Akash Prakasan2

All the cells in the body need to receive food (nutrients, metabolic products) and to dispose of waste products. The responsible system for that is cardiovascular system. It is responsible to supply food through the arteries and to return waste products through the veins for all living cells in the human body. This task is reached by a circulating fluid, the blood. The central location which all lines of supply originate from and return to is a small, very small, pump, the heart. The heart keeps the fluid in circulation. In the heart, there are two pumps, propelling blood into the pulmonary and systemic circulation and are combined into a single muscular organ to synchronously beat. Any disruption in the blood flow causes a disruption in food supply. Life is not possible without blood, but in the truth life is not possible without the circulation of blood. It must pump at all times, which it does by contracting and relaxing in a rhythmic pattern, approximately once every second, more than 86 thousand times every day, and about 2 billion times in a lifetime of 75 years, nonstop (Zamir, 2005). The blood ejected by the heart follows in the direction the arterial tree. Along the arterial tree, the arteries successively decrease in size, increase in number, undergo structural changes, and finish in arterioles that are as little as 10 m in diameter. The structure of the artery is quite complex. The main components of the vessel wall are endothelium, smooth muscle cells, elastic tissue, collagen, and connective tissue. The arteries are targets for diseases such as atherosclerosis or aneurysms that each year claims the lives of scores of people worldwide. The cardiovascular disease may be triggered or aggravated by mechanical stimuli, such as wall stress or stretch resulting from the blood pressure, or shear stress resulting from the blood flow (Wernig and Xu, 2002). Arteries can also adapt to long-term physiological conditions by thinning or thickening the muscular layer, and altering the relative composition and organization of the various assemblies of structural proteins in process generally know as remodelling. Bessa et al. (2011) showed that occurs remodelling in tail arterial bed from normotensive and hypertensive rats. As shown in Figure 1, the internal diameter of the proximal portion of the tail artery did not differ between Wistar rats and spontaneously hypertensive rats (SHR), whereas the diameter of the intermediate and distal portions of SHR tails arteries were significant

**1. Introduction** 

smaller than those of normotensive rats.

*1Department of Environmental Sciences and Technological* 

**Analysis in Aneurysm** 

*Rural Federal University of Semi-Arid 2Institute Dante Pazzanese of Cardiology* 

