# Computational Fluid Dynamics Analysis of Pulsatile Blood Flow Behavior in Modelled Stenosed Vessels with Different Severities

Authors: Mohsen Mehrabi; Saeed Setayeshi
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This study focuses on the behavior of blood flow in the stenosed vessels. Blood is modelled as an incompressible non-Newtonian fluid which is based on the power law viscosity model. A numerical technique based on the finite difference method is developed to simulate the blood flow taking into account the transient periodic behaviour of the blood flow in cardiac cycles. Also, pulsatile blood flow in the stenosed vessel is based on the Womersley model, and fluid flow in the lumen region is governed by the continuity equation and the Navier-Stokes equations. In this study, the stenosis shape is cosine by using Tu and Devil model. Comparing the results obtained from three stenosed vessels with 30%, 50%, and 75% area severity, we find that higher percent-area severity of stenosis leads to higher extrapressure jumps and higher blood speeds around the stenosis site. Also, we observe that the size of the stenosis in stenosed vessels does influence the blood flow. A little change on the cross-sectional value makes vast change on the blood flow rate. This simulation helps the people working in the field of physiological fluid dynamics as well as the medical practitioners.


---

## Body

## 1. Introduction

Cardiovascular disease is one of the major causes of deaths in developed countries. Most cases are associated with some form of abnormal flow of blood in stenotic arteries. Therefore, blood flow analysis through stenosed vessels has been identified as one of the important area of research in the recent few decades. One of the motivations to study the blood flow was to understand the conditions that may contribute to high blood pressure. Past studies indicated that one of the reasons a person having hypertension is when the blood vessel becomes narrow. Thus, many researches have been done for analysis of blood flow in stenosed vessels. Some recent investigations in these areas are cited in [1, 2]. Also, few investigations have been made to study the blood flow in arteries by assuming blood as Newtonian fluid [3, 4], but this assumption is not valid when there is a low shear rate. Taking into account the non-Newtonian behavior of blood, described by the micropolar fluid, Mekheimer and Kot [5] studied the micropolar fluid model for blood flow through a tapered artery with a stenosis. Analysis of nonlinear blood flow in a stenosed flexible artery was studied by [6]. They observed that arteries are narrowed by the development of atherosclerotic plaques that protrude into the lumen resulting in stenosed arteries. In addition, the pulsatile flow of blood through stenosed vessels is studied by Sankar and Lee [7]. Although blood flow has been modelled by many researchers, there have only been a few numerical studies on the flow in stenotic arteries using the realistic pulsatile flow conditions on the inlet and outlet [8]. Thus, the non-Newtonian fluid and unsteady flow properties for modeling of blood flow are two needed conditions that elimination of any of them prevents to obtain real results [9]. In this paper a simulation model is developed to study the unsteady state blood flow through a stenotic artery of different severity. Blood is modelled as a non-Newtonian fluid. Using the straight tube having three different sizes of stenosis, 30%, 50%, and 70%, numerical simulations are carried out for the flow field based on the finite difference method. Also, Marker and Cell (MAC) method [10] has been used to solve the governing equations of motion. In the end, dependence of the flow on the severity of stenosis is investigated.


## 2. Method

### 2.1. Modeling of Stenosed Vessel Geometry

The simulation model of the stenosed vessel is depicted in Figure1; a cosine-shaped vessel segment with an axially symmetric stenosis is modeled as a rigid tube with a circular cross-section. Let z-axis be taken along the axis of artery while r is the radial coordinates. The geometry of the stenosis is shown in Figure 2 and described as [11]
(2.1)r0(z)=a0-h2(1+cos(πzz0)),z∈[-z0,z0],
where r0(z) denotes the radius of cosine-shaped arterial segment in the constricted region. a0 is the constant radius of the straight artery in nonstenotic region. 2z0 is the axial length of the stenosis, and h is the measure of the degree of the stenosis.Figure 1
The geometric model of the stenosis.Figure 2
The 3D simulation model of the stenosed vessel.


### 2.2. Modeling of Pulsatile and Unsteady Blood Flow

In order to model the pulsatile blood flow, Womersley-Evans theory is used to obtain the spatial and temporal distribution of velocity profile [12]. The Womersley-Evans theory is defined as [12]
(2.2)ϑ(t,r0(z)a0)=2ϑ0(1-(r0(z)a0)2)+∑m=1∞|Vm||ψm(r0(z)a0,τm)|cos(mωt-∅m+χm(r0(z)a0,τm)),ϑm(t,r0(z)a0)=|Vm||ψm(r0(z)a0,τm)|cos(ωmt-∅m+χm),ψm(r0(z)a0,τm)=τmJ0(τm)-τmJ0(r0(z)a0τm)τmJ0(τm)-2J1(τm),χm(r0(z)a0,τm)=∠ψ(r0(z)a0,τm),τm=j3/2a0ρμωm,Womersleynumber=a0ρμωm,
where a0 is vessel radius, r0(z) is radial distance from vessel axis, ρ is blood density, μ is blood viscosity, m is number of harmonics in Fourier series function, and J0 and J1 are Bessel functions of zero and one degree. The simulation result of pulsatile blood flow is shown in Figure 3.Figure 3
Axial centerline velocity waveform obtained from Womersley model (atz=0).


### 2.3. Blood as a Non-Newtonian Fluid

Blood is non-Newtonian fluid, because of presence of various cells. This means that when shear stress is plotted against the shear rate at a given temperature, the plot shows a nonstraight line with a nonconstant slope as shown in Figure4. This slope is called the viscosity of the fluid. With the change of shear stress rate, the non-Newtonian fluid viscosity would be changed [13].Figure 4
Classification of non-Newtonian fluids.In this study, power law fluid model is used to simulate the behavior of blood fluid. A power law fluid, or the Ostwald-de Waele relationship, is a type of generalized Newtonian fluid for which the shear stress,τ, is given by velocity gradient of power “n” with μ being the flow consistency index, ∂u/∂y the shear rate or the velocity gradient perpendicular to the plane of shear, and n the flow behaviour index. Power law fluids can be subdivided into three different types of fluids based on the value of their flow behaviour index (see Table 1). The quantity μeff=K(∂u/∂y)n-1 represents an apparent or effective viscosity as a function of the shear rate.Table 1
Different states of flow behaviour index.
n
<1
1
>1
Type of fluid
Pseudoplastic
Newtonian fluid
Dilatant


### 2.4. Formulation of the Problem

#### 2.4.1. Equations of Motion

The viscous, incompressible flow in a long tube with stenosis at the specified position is considered. Let (r*,θ*,z*) be the cylindrical polar coordinates with z*-axis along the axis of symmetry of the tube. The region of interest is 0≤r*≤r0(z*), 0≤z*≤L* (L* being the finite length). Let u* and v* be the axial and the radial velocity components, respectively, p* the pressure, and ρ the density, and ν denotes the kinematic viscosity of the fluid. The fluid is assumed to be nonhomogeneous and incompressible, and the flow is also laminar. Blood in physiological conditions may be considered as incompressible [14]. We introduce the nondimensional variables t=t*U/a0, r=r*/a0, z=z*/a0, r0(z)=r0*(z*/a0)/a0, u=u*/U, ν=ν*/U, p=p*/ρU2,  where a0 is the radius of the straight portion of the tube and U is the maximum inlet speed of the fluid. The governing unsteady Navier-Stokes equations for incompressible fluid flow representing conservation of mass and momentum fluxes may be expressed in dimensionless variables as
(2.3)r∂u∂z+∂νr∂r=0,∂u∂t+∂uν∂r+∂u2∂z+uνr=-∂p∂z+1Re[cν(∂2u∂r2+1r∂u∂r+∂2u∂z2)+∂cν∂r(∂u∂r+∂ν∂z)],∂ν∂t+∂uν∂z+∂ν2∂r+ν2r=-∂p∂z+1Re[cν(∂2ν∂r2+1r∂ν∂r+∂2ν∂z2-νr2)+2∂cν∂r∂ν∂r],
where Re=Ua0/ν denotes the Reynolds number and cν=1+k[1-(r*/a0)n]=1+k[1-rn].


#### 2.4.2. Boundary Conditions

There is no shear along the axis of the tube which may be stated mathematically as(2.4)∂u∂r=0,ν=0,onr=0.
On the inner wall of the tube, the no-slip boundary conditions are
(2.5)u=ν=0,atr=r0(z).
At the inlet, the flow is considered to be fully developed and laminar, therefore:
(2.6)u=(1-r2),ν=0,atz=0.
While at the outlet, zero velocity gradient boundary conditions are imposed: ∂u/∂z=0, ∂ν/∂z=0 at z=L.


#### 2.4.3. Finite Difference Formulation

Finite-difference discretization of the equations has been carried out in the present work in staggered grid, popularly known as MAC cell. In this type of grid alignment, the velocities and the pressure are evaluated at different locations of the control volume. The time derivative terms are differenced according to the first order accurate two-level forward time differencing formula. The convective terms in the momentum equations are differenced with a hybrid formula consisting of central differencing and second order upwinding. The diffusive terms are differenced using the three-point central difference formula. The source terms are centrally differenced keeping the position of the respective fluxes at the centers of the control volumes. The pressure derivatives are represented by forward difference formulae. Discretization of the continuity equation at (i,j) cell delivers
(2.7)Rjr0(zi)ui,jn-ui-1,jnδz-Rj2∂r0(zi)∂zutc-ubcδR+Rljνi,jn-Rlj-1νi,j-1nδR=0,
where utc and ubc are given in Layek et al. [15].Herer0(zi) is calculated at z=zi, ∂r0(zi)/∂z denotes the derivative of r0(z) with respect to z and calculated at z=zi. The quantities (zi,Rj), (zli,Rlj) are the coordinates of the cell centre and the right top corner of the cell, respectively. Considering the source, convective, and diffusive terms at the nth time level, the momentum equation in z-direction in finite difference form may be put as
(2.8)ui,j(n+1)-ui,jnδt=pi,jn-pi+1,jnδz+Rjr0(zli)∂r0(zli)∂zpt-pbδR+Ucdi,jn,
where the terms pt, pb and Ucdi,jn are defined in Layek et al. [15] and pt, pb stand for pressure at the top and bottom middle positions of the u-momentum equation at nth time level at (i,j)th cell. The diffusive terms are discretized centrally, and central difference formula is used for the mixed derivative ∂2u/∂z∂R in uniform grid sizes.The finite difference equation approximating the momentum equation in theR-direction is
(2.9)νi,j(n+1)-νi,jnδt=1r0(z)pi,jn-pi,j+1nδR+Vcdi,jn,
where Vcdi,jn=(1/Re)Diffνi,jn-Conνi,jn.HereVcdi,jn is the discretization of convective and diffusive terms of ν-momentum equation at the nth time level at cell (i,j). The diffusive and the convective terms in the ν-momentum equation are differenced similar to that in u-momentum for the convective flux. The Poisson equation for pressure is obtained by combining the discretized form of the momentum and continuity equations. The final form of the Poisson equation for pressure is
(2.10)(A+B+C+D)pi,jn-Api+1,jn-Bpi-1,jn+A1pi+1,j+1n-A1pi+1,j-1n-A2pi-1,j+1n+A2pi-1,j-1n-(C-A1+A2)pi,j+1n-(D+A1-A2)pi,j-1n=-[Diνi,jnδt+Rjr0(zi)Ucdi,jn-Ucdi-1,jnδz+RljVcdi,jn-Rlj-1Vcdi,j-1nδR],
where A,B,C,D,A1,A2 are all given in Layek et al. [15]. Here Diνi,jn is the finite-difference representation of the divergence of the velocity field at cell (i,j).


## 2.1. Modeling of Stenosed Vessel Geometry

The simulation model of the stenosed vessel is depicted in Figure1; a cosine-shaped vessel segment with an axially symmetric stenosis is modeled as a rigid tube with a circular cross-section. Let z-axis be taken along the axis of artery while r is the radial coordinates. The geometry of the stenosis is shown in Figure 2 and described as [11]
(2.1)r0(z)=a0-h2(1+cos(πzz0)),z∈[-z0,z0],
where r0(z) denotes the radius of cosine-shaped arterial segment in the constricted region. a0 is the constant radius of the straight artery in nonstenotic region. 2z0 is the axial length of the stenosis, and h is the measure of the degree of the stenosis.Figure 1
The geometric model of the stenosis.Figure 2
The 3D simulation model of the stenosed vessel.


## 2.2. Modeling of Pulsatile and Unsteady Blood Flow

In order to model the pulsatile blood flow, Womersley-Evans theory is used to obtain the spatial and temporal distribution of velocity profile [12]. The Womersley-Evans theory is defined as [12]
(2.2)ϑ(t,r0(z)a0)=2ϑ0(1-(r0(z)a0)2)+∑m=1∞|Vm||ψm(r0(z)a0,τm)|cos(mωt-∅m+χm(r0(z)a0,τm)),ϑm(t,r0(z)a0)=|Vm||ψm(r0(z)a0,τm)|cos(ωmt-∅m+χm),ψm(r0(z)a0,τm)=τmJ0(τm)-τmJ0(r0(z)a0τm)τmJ0(τm)-2J1(τm),χm(r0(z)a0,τm)=∠ψ(r0(z)a0,τm),τm=j3/2a0ρμωm,Womersleynumber=a0ρμωm,
where a0 is vessel radius, r0(z) is radial distance from vessel axis, ρ is blood density, μ is blood viscosity, m is number of harmonics in Fourier series function, and J0 and J1 are Bessel functions of zero and one degree. The simulation result of pulsatile blood flow is shown in Figure 3.Figure 3
Axial centerline velocity waveform obtained from Womersley model (atz=0).


## 2.3. Blood as a Non-Newtonian Fluid

Blood is non-Newtonian fluid, because of presence of various cells. This means that when shear stress is plotted against the shear rate at a given temperature, the plot shows a nonstraight line with a nonconstant slope as shown in Figure4. This slope is called the viscosity of the fluid. With the change of shear stress rate, the non-Newtonian fluid viscosity would be changed [13].Figure 4
Classification of non-Newtonian fluids.In this study, power law fluid model is used to simulate the behavior of blood fluid. A power law fluid, or the Ostwald-de Waele relationship, is a type of generalized Newtonian fluid for which the shear stress,τ, is given by velocity gradient of power “n” with μ being the flow consistency index, ∂u/∂y the shear rate or the velocity gradient perpendicular to the plane of shear, and n the flow behaviour index. Power law fluids can be subdivided into three different types of fluids based on the value of their flow behaviour index (see Table 1). The quantity μeff=K(∂u/∂y)n-1 represents an apparent or effective viscosity as a function of the shear rate.Table 1
Different states of flow behaviour index.
n
<1
1
>1
Type of fluid
Pseudoplastic
Newtonian fluid
Dilatant


## 2.4. Formulation of the Problem

### 2.4.1. Equations of Motion

The viscous, incompressible flow in a long tube with stenosis at the specified position is considered. Let (r*,θ*,z*) be the cylindrical polar coordinates with z*-axis along the axis of symmetry of the tube. The region of interest is 0≤r*≤r0(z*), 0≤z*≤L* (L* being the finite length). Let u* and v* be the axial and the radial velocity components, respectively, p* the pressure, and ρ the density, and ν denotes the kinematic viscosity of the fluid. The fluid is assumed to be nonhomogeneous and incompressible, and the flow is also laminar. Blood in physiological conditions may be considered as incompressible [14]. We introduce the nondimensional variables t=t*U/a0, r=r*/a0, z=z*/a0, r0(z)=r0*(z*/a0)/a0, u=u*/U, ν=ν*/U, p=p*/ρU2,  where a0 is the radius of the straight portion of the tube and U is the maximum inlet speed of the fluid. The governing unsteady Navier-Stokes equations for incompressible fluid flow representing conservation of mass and momentum fluxes may be expressed in dimensionless variables as
(2.3)r∂u∂z+∂νr∂r=0,∂u∂t+∂uν∂r+∂u2∂z+uνr=-∂p∂z+1Re[cν(∂2u∂r2+1r∂u∂r+∂2u∂z2)+∂cν∂r(∂u∂r+∂ν∂z)],∂ν∂t+∂uν∂z+∂ν2∂r+ν2r=-∂p∂z+1Re[cν(∂2ν∂r2+1r∂ν∂r+∂2ν∂z2-νr2)+2∂cν∂r∂ν∂r],
where Re=Ua0/ν denotes the Reynolds number and cν=1+k[1-(r*/a0)n]=1+k[1-rn].


### 2.4.2. Boundary Conditions

There is no shear along the axis of the tube which may be stated mathematically as(2.4)∂u∂r=0,ν=0,onr=0.
On the inner wall of the tube, the no-slip boundary conditions are
(2.5)u=ν=0,atr=r0(z).
At the inlet, the flow is considered to be fully developed and laminar, therefore:
(2.6)u=(1-r2),ν=0,atz=0.
While at the outlet, zero velocity gradient boundary conditions are imposed: ∂u/∂z=0, ∂ν/∂z=0 at z=L.


### 2.4.3. Finite Difference Formulation

Finite-difference discretization of the equations has been carried out in the present work in staggered grid, popularly known as MAC cell. In this type of grid alignment, the velocities and the pressure are evaluated at different locations of the control volume. The time derivative terms are differenced according to the first order accurate two-level forward time differencing formula. The convective terms in the momentum equations are differenced with a hybrid formula consisting of central differencing and second order upwinding. The diffusive terms are differenced using the three-point central difference formula. The source terms are centrally differenced keeping the position of the respective fluxes at the centers of the control volumes. The pressure derivatives are represented by forward difference formulae. Discretization of the continuity equation at (i,j) cell delivers
(2.7)Rjr0(zi)ui,jn-ui-1,jnδz-Rj2∂r0(zi)∂zutc-ubcδR+Rljνi,jn-Rlj-1νi,j-1nδR=0,
where utc and ubc are given in Layek et al. [15].Herer0(zi) is calculated at z=zi, ∂r0(zi)/∂z denotes the derivative of r0(z) with respect to z and calculated at z=zi. The quantities (zi,Rj), (zli,Rlj) are the coordinates of the cell centre and the right top corner of the cell, respectively. Considering the source, convective, and diffusive terms at the nth time level, the momentum equation in z-direction in finite difference form may be put as
(2.8)ui,j(n+1)-ui,jnδt=pi,jn-pi+1,jnδz+Rjr0(zli)∂r0(zli)∂zpt-pbδR+Ucdi,jn,
where the terms pt, pb and Ucdi,jn are defined in Layek et al. [15] and pt, pb stand for pressure at the top and bottom middle positions of the u-momentum equation at nth time level at (i,j)th cell. The diffusive terms are discretized centrally, and central difference formula is used for the mixed derivative ∂2u/∂z∂R in uniform grid sizes.The finite difference equation approximating the momentum equation in theR-direction is
(2.9)νi,j(n+1)-νi,jnδt=1r0(z)pi,jn-pi,j+1nδR+Vcdi,jn,
where Vcdi,jn=(1/Re)Diffνi,jn-Conνi,jn.HereVcdi,jn is the discretization of convective and diffusive terms of ν-momentum equation at the nth time level at cell (i,j). The diffusive and the convective terms in the ν-momentum equation are differenced similar to that in u-momentum for the convective flux. The Poisson equation for pressure is obtained by combining the discretized form of the momentum and continuity equations. The final form of the Poisson equation for pressure is
(2.10)(A+B+C+D)pi,jn-Api+1,jn-Bpi-1,jn+A1pi+1,j+1n-A1pi+1,j-1n-A2pi-1,j+1n+A2pi-1,j-1n-(C-A1+A2)pi,j+1n-(D+A1-A2)pi,j-1n=-[Diνi,jnδt+Rjr0(zi)Ucdi,jn-Ucdi-1,jnδz+RljVcdi,jn-Rlj-1Vcdi,j-1nδR],
where A,B,C,D,A1,A2 are all given in Layek et al. [15]. Here Diνi,jn is the finite-difference representation of the divergence of the velocity field at cell (i,j).


## 2.4.1. Equations of Motion

The viscous, incompressible flow in a long tube with stenosis at the specified position is considered. Let (r*,θ*,z*) be the cylindrical polar coordinates with z*-axis along the axis of symmetry of the tube. The region of interest is 0≤r*≤r0(z*), 0≤z*≤L* (L* being the finite length). Let u* and v* be the axial and the radial velocity components, respectively, p* the pressure, and ρ the density, and ν denotes the kinematic viscosity of the fluid. The fluid is assumed to be nonhomogeneous and incompressible, and the flow is also laminar. Blood in physiological conditions may be considered as incompressible [14]. We introduce the nondimensional variables t=t*U/a0, r=r*/a0, z=z*/a0, r0(z)=r0*(z*/a0)/a0, u=u*/U, ν=ν*/U, p=p*/ρU2,  where a0 is the radius of the straight portion of the tube and U is the maximum inlet speed of the fluid. The governing unsteady Navier-Stokes equations for incompressible fluid flow representing conservation of mass and momentum fluxes may be expressed in dimensionless variables as
(2.3)r∂u∂z+∂νr∂r=0,∂u∂t+∂uν∂r+∂u2∂z+uνr=-∂p∂z+1Re[cν(∂2u∂r2+1r∂u∂r+∂2u∂z2)+∂cν∂r(∂u∂r+∂ν∂z)],∂ν∂t+∂uν∂z+∂ν2∂r+ν2r=-∂p∂z+1Re[cν(∂2ν∂r2+1r∂ν∂r+∂2ν∂z2-νr2)+2∂cν∂r∂ν∂r],
where Re=Ua0/ν denotes the Reynolds number and cν=1+k[1-(r*/a0)n]=1+k[1-rn].


## 2.4.2. Boundary Conditions

There is no shear along the axis of the tube which may be stated mathematically as(2.4)∂u∂r=0,ν=0,onr=0.
On the inner wall of the tube, the no-slip boundary conditions are
(2.5)u=ν=0,atr=r0(z).
At the inlet, the flow is considered to be fully developed and laminar, therefore:
(2.6)u=(1-r2),ν=0,atz=0.
While at the outlet, zero velocity gradient boundary conditions are imposed: ∂u/∂z=0, ∂ν/∂z=0 at z=L.


## 2.4.3. Finite Difference Formulation

Finite-difference discretization of the equations has been carried out in the present work in staggered grid, popularly known as MAC cell. In this type of grid alignment, the velocities and the pressure are evaluated at different locations of the control volume. The time derivative terms are differenced according to the first order accurate two-level forward time differencing formula. The convective terms in the momentum equations are differenced with a hybrid formula consisting of central differencing and second order upwinding. The diffusive terms are differenced using the three-point central difference formula. The source terms are centrally differenced keeping the position of the respective fluxes at the centers of the control volumes. The pressure derivatives are represented by forward difference formulae. Discretization of the continuity equation at (i,j) cell delivers
(2.7)Rjr0(zi)ui,jn-ui-1,jnδz-Rj2∂r0(zi)∂zutc-ubcδR+Rljνi,jn-Rlj-1νi,j-1nδR=0,
where utc and ubc are given in Layek et al. [15].Herer0(zi) is calculated at z=zi, ∂r0(zi)/∂z denotes the derivative of r0(z) with respect to z and calculated at z=zi. The quantities (zi,Rj), (zli,Rlj) are the coordinates of the cell centre and the right top corner of the cell, respectively. Considering the source, convective, and diffusive terms at the nth time level, the momentum equation in z-direction in finite difference form may be put as
(2.8)ui,j(n+1)-ui,jnδt=pi,jn-pi+1,jnδz+Rjr0(zli)∂r0(zli)∂zpt-pbδR+Ucdi,jn,
where the terms pt, pb and Ucdi,jn are defined in Layek et al. [15] and pt, pb stand for pressure at the top and bottom middle positions of the u-momentum equation at nth time level at (i,j)th cell. The diffusive terms are discretized centrally, and central difference formula is used for the mixed derivative ∂2u/∂z∂R in uniform grid sizes.The finite difference equation approximating the momentum equation in theR-direction is
(2.9)νi,j(n+1)-νi,jnδt=1r0(z)pi,jn-pi,j+1nδR+Vcdi,jn,
where Vcdi,jn=(1/Re)Diffνi,jn-Conνi,jn.HereVcdi,jn is the discretization of convective and diffusive terms of ν-momentum equation at the nth time level at cell (i,j). The diffusive and the convective terms in the ν-momentum equation are differenced similar to that in u-momentum for the convective flux. The Poisson equation for pressure is obtained by combining the discretized form of the momentum and continuity equations. The final form of the Poisson equation for pressure is
(2.10)(A+B+C+D)pi,jn-Api+1,jn-Bpi-1,jn+A1pi+1,j+1n-A1pi+1,j-1n-A2pi-1,j+1n+A2pi-1,j-1n-(C-A1+A2)pi,j+1n-(D+A1-A2)pi,j-1n=-[Diνi,jnδt+Rjr0(zi)Ucdi,jn-Ucdi-1,jnδz+RljVcdi,jn-Rlj-1Vcdi,j-1nδR],
where A,B,C,D,A1,A2 are all given in Layek et al. [15]. Here Diνi,jn is the finite-difference representation of the divergence of the velocity field at cell (i,j).


## 3. Results and Discussion

In this study we have simulated and analyzed pulsatile flow of a power law fluid as a model for blood flow in the cardiovascular system. The model is used to study the critical flow in stenotic arteries with three different severities of 30%, 50%, and 70% (shown in Figures8, 11, and 14). The solutions were computed for five cardiac cycles to ensure reproducibility of the pulsed characteristic flow. The time step size Δt used for each model is allowed to have different values. For the case with more percentage area severity, we use smaller time step size. For the case of 30% area severity, the minimum time step is taken to be Δtmin=0.005s, the maximum time step is Δtmax=0.01s, and for 50% and 70% area severity, the time steps are changed to Δtmin=0.001s, Δtmax=0.005s.Upstream from the stenosis, the velocity profile in thez-direction is parabolic as shown in Figure 5, and the fluid passes through the stenosis at high speed, especially at the throat of the stenosis. Downstream from the stenosis region, the distal part, the flow has stair-step shape profile, and the longitudinal velocity uz is negative (along the negative z direction) in the recirculation region. A region of reversal flow occurs at the downstream, next to the stenosis whereas the jet impinging occurs at the throat of the stenosis. Higher area blockage severity leads to larger pressure dropping around the stenosis and consequently gives higher speed in the stenosis area.The velocityuz at the peak of systole in the Oxy plane of the 70% stenotic tube: (a) at the upstream cross-section z=1.5 cm, (b) at the throat cross-section z=2.5 cm, and (c) at the downstream cross-section z=2.7 cm.
(a)
(b)
(c)A rapid fall in pressure is observed as the occlusion is approached. Higher percentage area severity leads to greater pressure drops around the stenosis (shown in Figures9, 12, and 15). It also leads to higher speeds in the stenosis area (shown in Figures 10, 13, and 16). The results clearly show the relation between pressure and velocity field. These figures illustrate the pressure distribution along a longitudinal line and the flow velocity of blood at peak systole related to stenosed vessel with different severities. It shows that the pressure drops very quickly near the stenosis site and creates a jet flow at the throat of the stenosis. The flow is accelerating when a negative pressure gradient exists at the stenosis site. The flow decelerates when an adverse pressure gradient exists. This extrapressure jump helps to impel the flow passing through the narrowing channel. Figure 6 shows the time characteristics of the pulsatile flow in terms of mean flow, pressure, and shear rate. In a healthy artery, the wall shear stress is approximately 15 dyn/cm2 [16]. To determine the critical flow condition, measurement of wall shear stress using numerical experiments becomes necessary. The shear rate and wall shear stress were computed for a model with 70% area severity. The magnitude of shear rate along a longitudinal line increases sharply before the occlusion is approached and reaches a maximum value near the center of the throat around 2×104s-1 at the peak of the systole and at 8×103s-1 at the peak of the diastole. It then decreases in the downstream. To depict the wall shear stress along the arterial wall, we plot the solution on the plane representing the wall surface where the stenosis is located at the center. The direction of the wall shear stress oscillates in the recirculation zone at downstream as shown in Figure 7.Pulsatile flow velocity, pulse pressure, and variation of shear rate with respect to time at an upstream point (2.3,0,0) for a 50% stenotic tube.
(a)
(b)
(c)Wall shear stress at the end of diastolet=3.70 s: (a) surface plot and (b) contour plot.
(a)
(b)Figure 8
Blood flow velocity distribution at peak systole obtained from simulation of stenosed vessel with severities of 30%.Figure 9
Pressure along a longitudinal line at peak systole related to stenosed vessel with severities of 30%.Figure 10
Flow velocity of blood along a longitudinal line at peak systole related to stenosed vessel with severities of 30%.Figure 11
Blood flow velocity distribution at peak systole obtained from simulation of stenosed vessel with severities of 50%.Figure 12
Pressure along longitudinal line at peak systole related to stenosed vessel with severities of 50%.Figure 13
Flow velocity of blood along a longitudinal line at peak systole related to stenosed vessel with severities of 50%.Figure 14
Blood flow velocity distribution at peak systole obtained from simulation of stenosed vessel with severities of 70%.Figure 15
Pressure along a longitudinal line at peak systole related to stenosed vessel with severities of 70%.Figure 16
Flow velocity of blood along a longitudinal line at peak systole related to stenosed vessel with severities of 70%.The results also show a similar pattern in the pulsatile velocity, in the pulse pressure, and in the variation of shear rate in cardiac cycles. These confirm the features of the characteristic of the periodic motion. Therefore, in the presence of a narrowing vessel lumen with different area severity, the flow experiences resistance, which causes an increase in the shear stress and in the pressure drop. Higher percent-area severity of stenosis produces a higher pressure drop, a higher blood speed, a higher shear rate, and a higher wall shear stress.Comparing the results obtained from three stenotic tubes with 30%, 50%, and 70% area severity, we find that higher percent area severity of stenosis leads to higher extrapressure jumps, higher blood speeds around the stenosis site, higher shear rate, and higher wall shear stress.


## 4. Conclusion

In this paper, we have derived a simple mathematical model that can represent the blood flow in the arteries. We observe that the size of the stenosis in stenosed vessels does influence the blood flow. A little change in the cross-sectional value makes vast change in the blood flow rate.It should be noted that blood flow in a small stenotic artery is an extremely complex phenomenon. There are many unresolved modeling problems such as the flow in the arterial wall which is deformed during the cardiac period. The presented work only focuses on blood flow in the stenosed vessels. Also, this simulation helps the people working in the field of physiological fluid dynamics as well as the medical practitioners.



---

*Source: 804765-2012-08-27.xml*