Pulsatile flow of non-Newtonian fluids through arterial stenoses

Pergamon 0021-9290(95)00151-4

J. Biomechanics, Vol. 29, No. 7, pp. 899-908, 1996 Copyright/~ 1996 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0021 9290/96 $15.00 + ,00

PULSATILE FLOW OF NON-NEWTONIAN FLUIDS THROUGH ARTERIAL STENOSES Cheng Tu* and Michel Devillet *Center for Systems Engineering and Applied Mechanics, Universit6 Catholique de Louvain B~timent Euler, avenue G. Lemaitre, 4 B-1348 Louvain-la-Neuve, Belgium; and ?IMHEF, Swiss Federal Institute of Technology, DGM-ECUBLENS, CH-I015 Lausanne, Switzerland Abstract--The problem of blood flow through stenoses is solved using the incompressible generalized Newtonian model. The Herschel-Bulkley, Bingham and power-law fluids are incorporated. The geometry corresponds to a rigid circular tube with a partial occlusion. Calculations are performed by a Galerkin finite-element method. For the pulsatile case, a predictor-corrector time marching scheme is used with an adaptive time step. Results are obtained for steady and pulsatile physiological flows. Computations show that the memory effects taken into account in the model affect deeply the flow compared with the Newtonian reference case. The disturbances are stronger by their vorticity intensity and persist after the geometrical obstacle. This is especially true for severe stenoses. Copyright ~) 1996 Elsevier Science Ltd. Keywords: Stenosis; Non-Newtonian fluid; Pulsatile flow.

INTRODUCTION There is considerable evidence that vascular fluid dynamics play an important role in the development and progression of arterial stenosis, one of the most widespread diseases in human beings. In order to understand the effect of stenosis on blood flow through and beyond the narrowed segment of the artery, many studies have been undertaken experimentally and theoretically. Until now, most theoretical analysis has been done assuming the b l o o d as a Newtonian fluid: (Cheng et al., 1973, 1974; Daly, 1976; Deshpande et al., 1976; Fukushima et al., 1982; Kawaguti and Hamano, 1983; Lee and Fung, 1970; Liou et al., 1981; Mc Donald, 1979; Morgan and Young, 1974; O'Brien and Ehrlich, 1985; Padmanabhan, 1980; Pollard, 1981; Tu et al., 1992; Tutty, 1992; Wille and Walloe, 1981). The assumption of Newtonian behaviour of blood is acceptable for high shear rate flow, e.g. in the case of flow through large arteries. It is not, however, valid when the shear rate is low (0.1 s -1) which is the case in small arteries and in the downstream side of the stenosis. The non-Newtonian character of blood is typical in small arteries and veins where the presence of the cells induces that specific behaviour. It has also been pointed out that in some diseased conditions, e.g. patients with severe myocardial infarction, cerebrovascular diseases and hypertension, blood exhibits remarkable non-Newtonian properties (Chien, 1981). Starting with Shukla et al. (1980), a number of researchers have studied flow of non-Newtonian fluids through arterial stenosis: Chaturani and Samy (1985, 1986), Srivastava (1985), Chakravarty (1987), Nakamura and Sawada (1988, 1990),

Received in final form 4 September 1995.

Address correspondence to: M. Deville.

Chakravarty and Datta (1989, 1992), and Misra et al. (1993). Nevertheless, to our knowledge, only three papers take into account both the unsteadiness and the nonNewtonian effects: Chaturani and Samy consider blood obeying the Casson's fluid model; Nakamura et al. use the biviscosity (modified Bingham) model; and Misra et aL apply a model in which blood is treated as a viscoelastic fluid with fluid motion equations derived on the basis of long-wave approximations. In these three papers, the velocity waveform is not a physiological one. Nakamura and Sawada use a rectangular waveform; Chaturani and Samy use a sinusoidal one; the analysis of Misra et al. is restricted to propagation of small-amplitude harmonic waves. Furthermore, the computation of Chaturani and Samy and Misra et al. is restricted to mild stenoses. The results presented in these papers show only the change of the separated region, the velocity distribution, the resistance and the wall shear stress; however, the influence of stenosis on pressure, streamlines and vorticity is not shown. Thus, there is much left to be studied with the simultaneous effects of physiological pulsatile flow and non-Newtonian behaviour of blood in the presence of stenosis. Faster computers and more advanced computational algorithms have now made it possible to perform studies using a finer mesh and under more realistic situations. In this study, the blood is treated using a HerschelBulkley model, which represents blood flow fairly closely. (Blood obeys Casson's model only for moderate shear rate flows 7 < 10 s -1, according to some studies, e.g. Fung, 1984.) The other advantage of the HerschelBulkley model is that it can be restricted to model the behaviour of Newtonian fluid, Bingham fluid and powerlaw fluid by taking the appropriate values of the parameters. We use a finite-element approximation for solving the unsteady generalized Navier-Stokes equations. The algorithm consists of a predictor-corrector scheme 899

900

C. Tu and M. Deville

in which the nonlinear system generated by the corrector step is solved using a one-step Newton method. It can automatically select the appropriate time step in such a way that the method is cost-efficient. We first study the effect of mild and severe stenoses on steady flow. The model gives results which agree with previous computations. We also perform simulations for a physiological pulsatile flow through a severe stenosis. The model is capable of predicting the hemodynamic features most interesting to the physiologists: the detailed pictures of the velocity profile, the formation of the vortex in the separated regions, the pressure drop across the stenosis, the wall shear stress as well as the vorticity contours. The rest of the paper is organized as follows. In the next section, we present the basic equations. This is followed by the method of solution. The last two sections are devoted, respectively, to results and conclusions.

BASIC EQUATIONS AND NUMERICAL METHOD

It is well known that blood behaves as a non-Newtonian fluid under certain flow conditions (see for example, Fung, 1984). Two characteristic features emerge from experimental data: (i) the presence of a yield stress and (ii) the dependence of the viscosity with respect to the shear rate. To take these effects into account, we use the Herschel-Bulkley law to model the fluid behaviour of blood flow (Scott Blair and Spanner, 1974; Kapur, 1985). The momentum equation is D¥

P~--7 = - Vp + V.T, L) t

(1)

and the incompressibility constraint imposes the divergence-free condition on the velocity v: V.v=O.

(2)

In equation (1), p is the density, p the pressure and T the extra stress tensor. The material time derivative is defined as usual by the relation D D'--t= t?t + v.V.

1

D = ~ (Vv + vvT),

q(33) =

{

Zy/~ -1-K33 " - 1

33-> 33c,

Zy(2-U%)/33o+K[(2-n)+(n-1)33/%]33~ -~

33<%, (6)

where Zy is the yield stress, n the power index, 33c the critical shear rate, K the consistency factor (Crochet et al., 1992). We note that equation (6) has special cases: n = 1 and zy = 0, represents the Newtonian case with K the Newtonian viscosity; n = 1 and zy # 0 produces the Bingham fluid; n # 1 and ry = 0 corresponds to the power-law fluid. Equations (1)-(6) are treated in cylindrical coordinates and for the sake of simplicity, it is assumed that the geometry is axisymmetric (~/~0 = 0) and that the azimuthal velocity component vanishes. Although the vascular wall is a viscoelastic material, it is considered as a rigid tube. The boundary conditions are the following. On the rigid no-slip wall, the velocity field vanishes: v=0.

(7)

On the symmetry axis, the normal component of the velocity and the shear stress vanish: v.n=0,

T,z = 0 .

(4) (5)

where the superscript T indicates the transpose. The shear rate 33 is obtained from the second invariant of the tensor D such that 33= (2D:D) 1/2. In equation (4), the non-Newtonian viscosity satisfies the following

(8)

In the previous relations, n denotes the unit outward normal vector, while r and z are the radial and axial space directions, respectively. At the inflow section, the radial velocity vanishes while the axial component is calculated on the basis of the flow rate and the type of fluid. At the outflow section, we assume that the fluid comes out of the computational domain at its own velocity. This is ensured by imposing natural boundary conditions or stress-free conditions: f . = f = 0,

(9)

where the indices n and t refer to normal (axial) and tangential (radial) components of force f defined as f=n.a,

(3)

The Herschel-Bulkley model is a generalized Newtonian fluid such that the extra stress T is linearly dependent on the rate of deformation tensor D. However, the viscosity r/is a function of the shear rate 33. We have T=2q(33)D,

relation:

(10)

and a is the full stress tensor given by a = - pI + T.

(11)

The set of equations (1)-(6) is nonlinear because of the advection in the inertial term and the power-law behaviour in the constitutive relationship. Consequently, one must resort to a numerical technique, which is the finiteelement method in our case. The computational domain is broken into quadrilateral elements where velocities and pressures are approximated by quadratic and linear Lagrangian interpolants, respectively. This discretization avoids the presence of spurious pressure modes and ensures an optimal spatial convergence of the method. A typical mesh around the stenosis is provided in Tu et al. (1992) Fig. 1.

Pulsatile flow of non-Newtonian fluids Inserting equations (4) and (5) in equation (1) and applying the Galerkin projection to a weak formulation of this new equation leads to a semidiscrete set of equations which can be cast in a matrix form:

10'

M/~ + A(x)x + F(x) = O,

10°

901

X,

(12)

where the global vector x includes velocity and pressure nodal values. The time derivative of x is denoted b y / c The matrix M is the mass matrix, while A (x) incorporates the stiffness matrix resulting from the generalized viscous term and the nonlinear advective operator. The time integration is performed by a predictor-corrector method with the Euler explicit forward scheme for the predictor and the Euler implicit backward formula for the correction leading to the following nonlinear algebraic system:

\.,,,., Xx,

"\ 10-1¸

10 "z

lif t

'

I

'

I

'

101

10 °

I

102

M(x n+l _ x n) = _ Atn+I[A(xn+I)x n+l + F(x"+I)],

(13) where the superscript indicates the time level of the solution and At,+ 1 = t,+ 1 - t,. It has been found that for generalized Newtonian fluids with a low power index, Newton's method does not converge in many circumstances. Thus, this equation is solved first by a few Picard iterations followed by Newton's method. The timemarching scheme allows a variable time step in order to accomodate the dynamics of the flow. The philosophy behind the algorithm consists in using the difference between predicted and corrected values in order to evaluate the next time step in such a way that the local time error is inferior to a prescribed a priori level. This algorithm was proposed by Gresho et al. (1980). The resulting algebraic system is solved by the frontal method, a direct method based on Gaussian elimination.

RESULTS

The basic equations and the numerical methods mentioned above are applied to both steady and pulsatile flows through an axisymmetric stenosis. Although our main target is to study effects of pulsatile flow, we feel that it is useful to present results for the steady cases as well. This is because most of the numerical simulations done in recent years were restricted to mild numerical stenosis, e.g. Shukla et al. (1980), Srivastava (1985), Chaturani and Samy, (1985, 1986). We would like to verify our model by comparing with their results and to see if their findings on mild stenosis can be generalized to severe stenosis. For steady flow, we present results for both mild (25 %) and severe (75%) stenosis. For the pulsatile flow, only flow through 75% stenosis is studied. The shape of the stenosis is given by

R BI4 29-7-0

1

otherwise,

Fig. 1. Viscosity profile with respect to the shear rate. Dashed line: Bingham fluid; dotted line: Herschel-Bulkley fluid; continuous line: power-law fluid.

where h is the constricted radius (radius in the narrowest part, the neck); R is the unconstricted radius. The origin z = 0 corresponds to the neck of the stenosis. The finiteelement grid consists of 750 elements with 6262 velocity and 816 pressure unknowns. Blood is treated as the Herschel-Bulkley model with viscosity # = 4 x 10- 3 Pa s and density p = 1060 kg m - 3. We perform simulations for the following four cases. Case 1 (n = 1, r r = 0) is a Newtonian fluid; Case 2 (n = 1, zr = 0.2) is a Bingham fluid; Case 3 (n = 0.9, Zy = 0.2) is a Herschel-Bulkley fluid; and Case 4 (n = 0.8, zr = 0) is a power-law fluid. Figure 1 shows the viscosity versus the shear rate. The dashed line corresponds to the Bingham fluid, the dotted line to the Herschel-Bulkley fluid and the continuous line to the power-law fluid. Steady flow We first perform the calculation on 25% and 75% stenoses for the four cases above. The Reynolds number based on the radius, the inlet average velocity and the Newtonian kinematic viscosity v is defined by the relation Re = 2RU/v. For the Newtonian case, it is 67. The results are represented by the following nondimensional quantities: v v* = - - , Vm

p* =

p 2R(Op/Sz)o

,

~w z*

= --,

Zwo

where Vmis the mean axial velocity, (Op/Oz)o, the pressure gradient, Zwo, the wall shear stress, all taken far away from the stenosis. The same values of Vm (hence same flow-rate), (@/~z)o and Zwo are used for both Newtonian and non-Newtonian flows. The axial velocity profiles at different axial locations are illustrated in Fig. 2. First note that results of Cases 1 and 2 are almost identical. This shows that a small

902

C. Tu and M. Deville 2.5

i

i

i

2.5

!

i

i

i

i

i

i

i

2.5

i

I

I

I

'case I'

-'case 2' .....

2 ................ ~ ' : ~ '

cas

e

2

1.52 ~~,~I'2~'LL;~':'~:-:.:.:.:.:.:.:.:<. .

3 .......

%1.5

1

1

%

~?,

'case 3' ...... 'case 4 ............

0.5

0.5

2 ......

'case

I

'case i' - 'case'C~aS4'~i ese .....

1.5 1

-~2, \~

0.5 I

I

I

I

0.2

0.4

0.6

0.8

0

0

I

l

0

J

l

0.2

l

l

l

0.4

l

~

0.6

0.

0

l

I

I

I

0.2

0.4

0.6

0.8

r/R (b) 6

I

!

I

I

I

I

I

I

4.5

I

4

% 3

'case

2

'case 3

2

' case

I

0

I

I

0.i

4

I

..... ...... ..........

I

l

0.2

0.3

,

,

.'case

3 %

3.5

,

I

I

......

.............. i i<

2

' ...........

1.5 1 0.5

-0.5'

I

I

.........."...,.. .... 'case'Case 2'i'..... -2.5 ",.. 'case 3' ..... %

1.5 1 0.5 0

0.4

3'

I

3

2.52

\"? ~,,

l

,

3.54 [~'['~'['~i'i'i~""" 'caSe, case 21'...... --

,c~iiii\

5

r/R

I

i

I

I

0.4 0.6 r/R

0.2

.5

r/R (d)

0

0.8

0.2

0.4 0.6 r/R

(,)

0.8

(f)

Fig. 2. Axial velocity profiles at different axial locations for 25% and 75% stenoses. (a) Far away from the stenosis; (b) 25% stenosis at z/D = 0; (c) 25% stenosis at z/D = 2; (d) 75% stenosis at z/D = 0; (e) 75% stenosis at z/D = 2; (f) 75% stenosis at z/D = 5.

%

2.4 2.35 2.3 2.25 2.2 2.15 2.1 2.05 2

!

I

I

6

I

' c a s e I' ' c a s e 2' ..... ' c a s e 3' ...... ~case 4' ...........

....................

'

~

I

I

cane

case case ~iase

4

I

1'

- -

2' ..... 3 ....... 4' ...........

3.5 3 2.5

.....

1.95 1.9 1 85 ........i -i0 -5

%

i

5.5 5 4.5

2

i

,

0

5 Z/D

...... i............ i

i0

15

:.:.:..-.:.r.

1.5

........

20

(a)

-I0

i

i

I

I

I

-5

0

5

I0

15

20

Z/D

(b)

Fig. 3. Axial velocity variation along the symmetric axis: (a) 25% stenosis; (b) 75% stenosis. change of the yield shear stress would not have much influence on the flow behaviour. Figure 2(a) shows that at the same flow rate, the velocity near the axis for the Newtonian case is slightly higher than Cases 3 and 4, while this relative behaviour is reversed near the wall. From Fig. 2(b) and (c) we see that this is also true at the other axial locations for the 25% stenosis. This agrees with the previous results obtained by Chaturani and Samy. At z/D = 2, the effects of the stenosis on the flow are almost imperceptible for all the cases. For the 75% stenosis, at z/D = 0 [Fig. 2(d)], the velocity for the Newtonian case is also higher than those of Cases 3 and 4. However, at z/D = 2 and 5 [Fig. 2(e) and (f)], this is no longer true. At z/D = 5, velocities of Cases 1 and 2 have attained the same values as at the inlet, while for the other two cases the velocities are still much higher than their values at the inlet. From the axial velocity variation along the symmetric axis shown in Fig. 3, we see clearly

that the Newtonian case has the highest velocity at all axial locations for the 25% stenosis. For 75% stenosis, the Newtonian case has also the highest peak, but for a range of the post-stenotic location, the non-Newtonian cases have higher velocities than the Newtonian case. It shows that the non-Newtonian stenotic flow recovers to its undisturbed state less quickly than the Newtonian flow, e.g. the velocity of Case 4 for 75% stenosis only gets back to its entry value at z/D = 13, while for Case 1, the velocity assumes its entry value at z/D = 3. Figure 4 shows the streamlines for four cases of 75% stenosis. Note that in this figure, only the part of the artery near the stenosis is shown in order to yield more detail. We observe the separation phenomenon in all cases. It occurs downstream of the throat section and reattaches at some distal location. The Newtonian case has the shortest recirculation zone. As the power index n decreases, the separation point moves slightly upstream

Pulsatile flow of non-Newtonian fluids

903

Case 1

Case 3

Case 4

Fig. 4. Streamlines for four cases of 75% stenosis.

30

i

25

" ~ i ~i

!

I

I

i

'case 'case 'c a s e

20

I' 2 ...... 3' i.].~.[.].i

.......

50 45 40

' ~ i ~ i

.35

T

~i

20253015

' ' ' ' c a s e i' ' c a s e 2 ...... 'case 3 . . . . . . .

~.'ii/" " ' ~

e 4 ............

10 5 0

10

-5 -I0

I

I

I

I

I

-5

0

5 Z/D

i0

15

20

]. .........

0

-

-5

i -5

-i0

i

.

i

i

i

i

0

5 Z/D

I0

15

(a)

20

(b)

Fig. 5. Pressure variation along the symmetric axis: (a) 25% stenosis; (b) 75% stenosis.

2.5

I

I

'case 'case 'case 'case

2 1.5

i

i' 2' 3' 4'

-..... ...... ...........

2O 15

-

~

'case

0

5 Z/D

3

......

i0

15

_

1 0.5 f............................................. 0 -0.5 -i(

!

I

l

I

I

-5

0

5 Z/D

i0

15

20

-I0

-5

20

(b)

Ca)

Fig. 6. Shear stress on the wall: (a) 25% stenosis; (b) 75% stenosis.

towards the throat and the length of the recirculation zone increases. (The streamlines for four cases of 25% stenosis show no separation phenomenon and are not presented here.) Figure 5 illustrates the pressure distributions at the centerline of the artery. Note that zero pressure at the outlet of the artery is specified because the pressure is determined only up to a constant in the model. In all cases, a pressure drop caused by the stenosis is observed. The pressure reaches a slight local minimum at the throat (represented by the kink in the pressure curve) and then recovers downstream from there. The Newtonian flow has the highest pressure drop, which agrees with the results of N a k a m u r a and Sawada (1980). F o r all cases of 25% stenosis, the pressure drop is much smaller than for the cases of 75% stenosis.

30

,

!

0.2

0.4

,

,

0.6

0.8

2O 15 10 5 0 -5 0

t

Fig. 7. Physiological velocity profile (Siouffi et al., 1984) used in the pulsatile simulation. The velocity v is the spatial average axial velocity with units expressed in cm s- 1.

904

C. T u a n d M. Deville

Wall shear stress, another quantity of considerable physiological interest, is shown in Fig. 6. In all the cases, the wall shear stress increases sharply reaching a peak value near the throat but slightly upstream. Downstream from the point of the maximum the wall shear stress decreases rapidly and reverses direction. The normalized peak values for 25% stenosis are 2.24, 2.25, 1.43, and 1.06 for Cases 1, 2, 3 and 4, respectively. For 75% stenosis, they are 24.42, 24.46, 17.53, 14.82 for Cases 1, 2, 3 and 4, respectively. These values for the

t=O.1

cases of 25% stenosis are much smaller than those of 75% stenosis. They decrease as n decreases, as observed by Shukla et al. (1980), and Chaturani and Samy (1985). Pulsatile flow

We now turn to our main results. In order to mimic physiological flow conditions, we used a pulse velocity with the time-varying function similar to the one used by Siouffi et al. (1984). The period of the heart beat is

~-.

t=0.2 t=0.3 t=0.4 t=0.5

~'-.~

t=0.6 t--0.7 t=0.8 t=0.9 t----1.0

(a) t=O.1 t=0.2 t:O.3

t=0.4 t----0.5

~.

~

~

~

_

__

t=0.6 t=0.7 t=0.8

t=0.9 t=l.0

(b) Fig. 8. I n s t a n t a n e o u s streamlines d u r i n g a cycle of pulsatile flow at intervals of 0.1 s: (a) Case 1; (b) Case 3.

Pulsatile flow of non-Newtonian fluids

minimum at t = 0.7 s, a vortex develops upstream of the stenosis. As the flow accelerates again during diastole, this upstream vortex disappears and a down stream vortex forms. The vortex increases its size when the flow reaches its second peak and becomes smaller when the flow decelerates during diastole. The strength of the vortex at the beginning of the accelerating phase (t = 0.1 s), at the peak (t = 0.3 s) during systole, and when the flow reverses direction (t = 0.7 s) is further shown in Fig. 9. Stronger vortices are seen for the non-Newtonian case at all three different times, especially when the inflow velocity is high (at t = 0.3 s). The behaviour of the axial velocity profile at the peak of the inflow (t = 0.3 s) is illustrated in Fig. 10. The shapes of the velocity profiles in the part of the artery near the stenosis for Cases 1 and 3 are shown by the velocity vectors in Fig. 10(a) and (b). Comparing the velocity profiles downstream from the stenosis shows that the Herschel-Bulkley model has stronger effects

supposed to be 1 s. This is shown in Fig. 7. The corresponding peak Reynolds number where the characteristic velocity is the maximum spatial average velocity over the heart beat period is 132.5. The Womersley number for the Newtonian case defined by ~ = Rx/~, where o9 is the angular frequency of the pulsatile motion, is 1.25. For the pulsatile flow simulation, we concentrate on the more interesting case, the 75% stenosis. Figure 8 shows the instantaneous streamlines during a cycle of pulsatile flow at intervals of 0.1 s for Cases 1 and 3. As in the cases of steady flow, only the part of the artery near the stenosis is shown. The flow behaviours in both cases are similar except that the non-Newtonian flow shows a more disturbed pattern. In both cases, the length of the recirculation zone increases as the flow accelerates during systole and decreases as the flow decelerates. At t = 0.6 s, no streamlines are shown because the values are too small corresponding to a nearly zero inflow velocity. When the inflow reaches its negative

f

905

~

(a)

m,

(b) Fig. 9. Vorticity contours at t = 0.1 s, t = 0.3 s and t = 0.7 s: (a) Case 1; (b) Case 3,

IW W j ~ ' ~ , W ' y

y

•

Y

•

•

•

•

r

F

(a)

(b)

5 4.5

,

,

,

,

i

,

,

4.5

,

3.5

4

3.5 3 2.5 2 1.5 1 0.5 0

3

I' - 2 ' ..... 3 ...... 4 " ..........

'case

'case 'case

'case

2.5 2 1.5 1

i

4

i

3.5

"-'r:':'::'2~ ....... ' c a s e 1' -~ , ~ " , , " ' . , ' c a s e 2' ..... ~'.%*',,~'.,.case 3 ....... ~ase 4' ..........

3 2.5

%

0.I

,

,

0.2

,

,

0.3 r/R

(c)

~

,

0.4

I

0.5

i

<.::,

2

i

i' - 2' ..... 3 ' .....

, ...........

1 0.5 0

0 ,

i

"''"'...... ' c a s e -'"--, "..,' c a s e "'-,, " ' c a s e

1.5

0.5 ,

i

-0.5

-0.5

0.2

0.4

0.6 r/R

(a)

0.8

1

I

I

I

!

0.2

0.4

0.6

0.8

r/R

(e)

Fig. 10. Axial velocity profiles at t = 0.3 s: (a) general picture of Case 1; (b) general picture of Case 3; (c) comparison of 4 cases at z/D = 0; (d) comparison of 4 cases at z/D = 5; (e) comparison of 4 cases at z/D = 15.

906

C. Tu and M. Deville ing that the profile is still strongly m o d i f i e d by the

than the N e w t o n i a n m o d e l due to shear-thinning effect: at the point where the N e w t o n i a n profile returns to its

stenosis. From the axial velocity distribution at the centre of the artery at three different times of the cardiac cycle shown in Fig. 11, we notice that the velocity assumes its entry level slower than in the steady case [compare with Fig. 3(b)]. Even when the peak value is much smaller than the steady case at t = 0.1 s, it takes longer for the flow to recover from its disturbed states. For Case 4, the velocity at z/D = 18 is still far from its

shape before the stenosis, the Herschel-Bulkley profile is still highly disturbed. The details of the profiles for the four cases at z/D=O, 5 and 15 are shown in Fig. 10(c)-(e). From these figures we see that the velocity profile is parabolic at the upstream section and it becomes flattened at the throat of the stenosis. D o w n s t r e a m from the stenosis, the velocity presents a

reversal zone near the wall which disappears further down from the stenosis. At the throat of the stenosis, the velocities of the four cases are quite similar. For

entry level.

a wide range of the axial location d o w n s t r e a m from the stenosis, the n o n - N e w t o n i a n cases have higher

O n l y the results for t = 0.1, 0.3 a n d 0.7 s are s h o w n here. D u r i n g the accelerating phase of systole, the pressure

velocities than the Newtonian case. At z/D = 5, the reversal zone near the wall is still seen for Cases 3 and 4. This zone persists as far as z/D = 15 for Case 4, indicat-

falls at the throat of the stenosis and recovers partially downstream of the stenosis, as in the case of the steady flow. The drop of the pressure is larger in periods of

1.8

i

1.6

L "case 'case •% 'case "~,'case

1.4

~1.2

i

5

i

I' - 2' ..... 3 ....... 4" ...........

The general nature of the pressure field is suggested by its variation along the symmetry axis shown in Fig. 12.

i

i

i

-CI.3

i

~ ' c a s e i' - ~--~^,qa s e 2 ......

4.5

~i~,..a.,

4

.

.

.

.

-C).4 -[). 5

.

-C). 6

%3.5 3

-C). 7

0.8

2.5

-C). 8

0.6

2

1

0.4 -I0

I

I

I

5 Z/D

I0

15

I

-5

0

-..... ...... ..........

-0.9 I

1.5 -i0

20

I' 'case 2' ' c a s e 3' ' c a s e 4'

'case

~ I

I

-5

0

5 Z/D

I

I

I0

15

-i

20

-I0

I

I

I

I

-5

5 Z/D

I0

15

{b)

{a)

20

(c)

Fig. 11. Axial velocity variation along the symmetric axis: (a) at t = 0.1 s; (b) at t = 0.3 s and (c) at t = 0.7 s.

14

l

12

~,-, "-~::

i0 ,

,

i

'case 'case

I!

'case

i

60

i

i

50 40 !

I' - 2 ...... 3 ......

illlliiii::: 4............ k

~ -..

~

~i~[ "" i

30~

i

1

I' - 2 ...... 3 ....... 4' ...........

0

i

'case 'case 'case

'case

...

-i0

! -5

i 0

i 5 Z/D

i I0

i 15

I

-2

-4 -10

-2

I

20 -3

,, ........

I

-i

i0 |

I

-5

.,

""*

-20 -i0

20

-5

I

I

I

5 Z/D

i0

15

20

-6 -10

.....

4]'case i' - I ' c a s e 2 ...... I ' c a s e 3' ...... / ' % 4 ' c a s e 4' ......... i

i

i

!

i

-5

0

5 Z/D

i0

15

(b)

20

(c)

Fig. 12. Pressure variation along the symmetric axis: (a) at t = 0.1 s; (b) at t = 0.3 s and (c) at t = 0.7 s.

I0

9( 8( 7( 6( 5( "J 4( 3( 2( i(

8

6

~4 2

l

'case 'case

-I0

-5

0

5 Z/D

(a)

I0

15

20

1

i

-..... ...... ...........

--.~

-I( -2(

-2

i' 2' 3' 4'

•c a s e

(

0

i

'case

-10

! -5

I

I

I

5 Z/D

i0

15

20

i

0.5

0

....... ~".~"

-0.5 -I .~-1.5 -2.5 -3 -3.5 -4 -4.5 -10

•case

'case "case 'case

-5

(b)

Fig. 13. Shear stress on the wall: (a) at t = 0.1 s; (b) at t = 0.3 s and (c) at t = 0.7 s.

i' - 2' ..... 3 ....... 4' ...........

|

I

i

5 Z/D

I0

15

(e)

20

Pulsatile flow of non-Newtonian fluids relatively larger flow (t = 0.3 s). It is negative when the inflow is in the reversed direction (t = 0.7 s). At these times the pressure is smaller for Cases 3 and 4 than for Cases 1 and 2. The wall shear stress developed in the pulsatile flow is presented in Fig. 13. Results for t = 0.1, 0.3 and 0.7 s are shown, as in Fig. 12. During the accelerating phase of systole, it rises sharply and then decreases quickly reversing its direction at the vicinity of the stenosis, as in the case of steady flow. It also reaches a maximum at the times of maximum flow rate. The normalized peak values at t = 0.1 s are 9.31, 9.88, 7.48 and 6.10 for Cases 1, 2, 3 and 4, respectively. At t = 0.3 s, they are 78.80, 80.11, 69.61 and 67.12 for Cases 1, 2, 3 and 4, respectively. When the flow reverses direction at t = 0.7 s, the picture of the stress is also reversed. It is negative except when it rises near the throat of the stenosis. It has a sharp negative peak immediately after the throat, with normalized peak values - 4 . 0 1 , - 3 . 9 8 , - 2 . 9 5 and - 2 . 3 8 for Cases 1, 2, 3 and 4, respectively. The above results cannot be compared directly with the previous numerical results because the only available numerical simulations on unsteady flow through stenoses use different velocity waveforms (see Introduction). Furthermore, the computation of Chaturani and Samy (1986) is restricted to mild stenosis and done using the Casson's fluid model. N a k a m u r a and Sawada (1990) use the biviscosity fluid model. Finally, Misra et al. (1993) treat blood as a viscoelastic fluid with fluid equations derived on the basis of long-wave approximations (Misra et al., 1989), and their computation is also restricted to a mild stenosis.

CONCLUSION It is observed from this study that the rheological properties of blood can significantly affect the flow phenomena. For the steady state case, the results indicate that the n o n - N e w t o n i a n effects weaken the distortion of flow pattern, pressure distribution and wall shear stress associated with the mild (25%) stenosis. As the height of the stenosis increases (75% stenosis), the velocity at the neck, the pressure drop and wall shear stress increases but this increase is less due to non-Newtonian behaviour of the blood. However, we notice that for n o n - N e w t o n i a n flow through 75% stenosis, the influence of the geometrical disturbance affects the flow over a longer axial range. In the pulsatile case, the Herschel-Bulkley model gives rise to stronger vorticity production at the throat, and the streamlines and velocity profiles indicate that the flow patterns remain in a disturbed state over a long distance. The pressure drop and the wall shear stress are smaller for the non-Newtonian cases, as in the steady flow. The results demonstrate that the model is capable of predicting the hemodynamic features most interesting to physiologists. It can be used to predict fast stenotic flow

907

patterns on an individual basis. It can also be used for studying other parametric effects. REFERENCES

Chakravarty, S. (1987) Effects of stenosis on the flow-behaviour of blood in an artery. Int. J. Engng Sci. 25, 1003-1016. Chakravarty, S. and Datta, A. (1989) Effects of stenosis on arterial rheology through a mathematical model. Mathl. Cornput. Modelling 12, 1601-1612. Chakravarty, S. and Datta, A. (1992) Dynamic response of stenotic blood flow in vivo. Mathl. Comput. Modelling 16, 3-20. Chaturani, P. and Samy R. P. (1985) A study of non-Newtonian aspects of blood flow through stenosed arteries and its applications in arterial diseases. Biorheology 22, 521-531. Chaturani, P. and Samy R. P. (1986) Pulsatile flow of Casson's fluid through stenosed arteries with applications to blood flow. Biorheology 23, 499-511. Cheng, L. C., Robertson, J. M. and Clark, M. E. (1973) Numerical calculations of plane oscillatory nonuniform flow II. Parametric study of pressure gradient and frequency with square wall obstacles. J. Biomechanics 6, 521-538. Cheng, L. C., Robertson, J. M. and Clark, M. E. (1974) Calculations of plane pulsatile flow past wall obstacles. Computers and Fluids 2, 363-380. Chien, S. (1981) Hemorheology in clinical medicine. Recent Advances in Cardiovascular Diseases 2 (Suppl.), 21-26. Crochet, M. J., Debbaut, B., Keunings, R. and Marchal, J. M. (1992) Polyflow: a multi-purpose finite-element program for continuous polymer flows. In Computer Modelling of Extrusion and Other Continuous Processes (Edited by Brien, K. O), pp. 25-50. Hanser, Munich. Daly, B. J. (1976) A numerical study of pulsatile flow through stenosed canine femoral arteries. J. Biomechanics 9, 465-475. Deshpande, M. D., Giddens, D. P. and Mabon, R. F. (1976) Steady laminar flow through modelled vascular stenoses. J. Biomechanics 9, 165-174. Fukushima, T., Azuma, T. and Matsuzawa, T. (1982) Numerical analysis of blood flow in the vertebral artery. J. biomech. Engng 104, 143-147. Fung, Y. C. (1984) Biodynamics: Circulation. Springer, New York. Gresho, P. M., Lee, R. L. and Sani, R. L. (1980) On the time dependent solution of the incompressible Navier-Stokes equations in two and three dimensions. In Recent Advances in Numerical Methods in Fluids (Edited by Taylor, C. T. and Morgan, M.), Vol. 1, pp. 27-80. Pineridge Press, Swansea, U.K. Kapur, J. N., (1985) Mathematical Models in Biology and Medicine. Affiliated East-West Press, India. Kawaguti, M. and Hamano, A. (1983) Numerical study on post-stenotic dilatation. Biorheology 20, 507-518. Lee, J. and Fung, Y. (1970) Flow in locally constricted tubes at low Reynolds numbers. J. appl. Mech. 37, 9-16. Liou, R. J., Clark, M. E., Robertson, J. M. and Cheng, L. C. (1981) Three-dimensional simulation of steady flow past a partial stenosis. J. Biomechanics 14, 325-337. McDonald, D. A. (1979) On steady flow through modelled vascular stenoses. J. Biomechanics 12, 13-20. Misra, J. C., Patra, M. K. and Misra S. C., (1993) A nonNewtonian fluid model for blood flow through arteries under stenotic conditions. J. Biomechanics 26, 1129-1141. Morgan, B. E. and Young, D. F. (1974) An integral method for the analysis of flow in arterial stenoses. Bull. Math. Biol. 36, 39-53. Nakamura, M. and Sawada, T. (1988) Numerical study on the flow of a non-Newtonian fluid through an axisymmetric stenosis. J. hiomech. Engng 110, 137 143. Nakamura, M. and Sawada, T. (1990) Numerical study on the unsteady flow of non-Newtonian fluid. J. biomech. Engng 112, 100-103.

908

C. Tu and M. Deville

O'Brien, V. and Ehrlich, L. W. I. (1985) Simple pulsatile flow in an artery with a constriction. J. Biomechanics 18, 117-127. Padmanabhan, N. (1980) Mathematical model of arterial stenosis. Med. Biol. Enon 9 Comput. 18, 281-286. Pollard, A. (1981) A contribution on the effects of inlet conditions when modelling stenoses using sudden expansions. J. Biomechanics 14, 349-355. Scott Blair, G. W. and Spanner, D. C. (1974) An Introduction to Biorheolooy. Elsevier, Amsterdam. Shukla, J. B., Parihar, R. S. and Rao, B. R. P. (1980) Effects of stenosis on non-Newtonian flow of the blood in an artery. Bull. Math. Biol. 42, 283-294.

Siouffi, M. Pelissier, R., Farahifar, D. and Rieu, R. (1984) The effect of unsteadiness on the flow through stenoses and bifurcations. J. Biomechanics 17, 299 315. Srivastava, L. M. (1985) Flow of couple stress fluid through stenotic blood vessels. J. Biomechanics 18, 479 485. Tu, C., Deville, M., Dheur, L. and Vanderschuren, L. (1992) Finite-element simulation of pulsatile flow through arterial stenosis. J. Biomechanics 25, 1141-1152. Tutty, O. R. (1992) Pulsatile flow in a constricted channel. J. Biomech. Eng. 114, 50-54. Wille, S. O. and Walloe, L, (1981) Pulsatile pressure and flow in arterial stenosis simulated in a mathematical model. J. biomed. Engn9 3, 17-24.