Biosciences Mathematical 190 (2004) 203220
www.elsevier.com/locate/mbs
Modeling of arterial stenosis and its applications to blood diseases
R.N. Pralhad
a
a,*
, D.H. Schultz
b
Faculty of Applied Mathematics, Institute of Armament Technology, Girinagar, Pune 411 025, Maharashtra, India b Department of Mathematical Sciences, University of Wisconsin, Milwaukee, WI 53201, USA Received 30 January 2001; received in revised form 11 November 2002; accepted 9 January 2004 Available online 20 May 2004
Abstract Blood flow in a stenosed tube has been modeled in the present studies. Blood flow is assumed to be represented by a couple stress fluid. Flow parameters such as velocity, resistance to flow, and shear stress distribution have been computed for different suspension concentrations (haematocrit), and for the blood diseases; polycythemia, plasma cell dyscrasias, and for Hb SS (sickle cell). The results have been compared with the case of normal blood and for other theoretical models. The importance of size effects in blood flow studies has been highlighted. 2004 Published by Elsevier Inc.
Keywords: Stenosis; Blood flow modeling; Viscosity; Cardiovascular disease; Couple stress fluids; Hemorheology
1. Introduction One of the leading causes of the deaths in the world is due to heart diseases, and the most commonly heard names among the same are ischemia, atherosclerosis, and angina pectoris. Ischemia is the deficiency of the oxygen in a part of the body, usually temporary. It can be due to a constriction (stenosis) or obstruction in the blood vessel supplying that part. Atherosclerosis is a type of arteriosclerosis. It comes from the Greek words athero (meaning gruel or paste) and sclerosis (hardness). It involves deposits of fatty substances, cholesterol, cellular waste products, calcium and fibrin (a clothing material in the blood) in the inner lining of an artery. The build up
*
Corresponding author. E-mail address: [e-mail] (R.N. Pralhad).
0025-5564/$ - see front matter 2004 Published by Elsevier Inc. doi:10.1016/j.mbs.2004.01.009
204
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
that results is called plaque. Plaque may partially or totally block the blood flow through an artery. Two things that can happen where plaque occurs are (i) bleeding (hemorrhage) into the plaque, and (ii) formation of a blood clot (thrombus) on the plaque's surface. If either of these occurs and blocks the entire artery, a heart attack or stroke may result. Atherosclerosis affects large and medium sized arteries. The type of artery where the plaque develops varies with each person. A symptom complex of ischemic heart disease characterized by paroxysmal attacks of chest pain, usually substernal or pre-cordial is referred as angina pectoris. Usually high-grade stenosis with acute coronary changes results in sudden cardiac arrest (or death) which strikes 300 000400 000 persons annually around the globe. Owing to its serious concern, a major research work is being done over all parts of the world for early detection and prevention from being affected by the cardiac attack and the development of the art therapies for the diagnosis of the heart diseases. In cardiac-related problems, the affected arteries get hardened as a result of accumulation of fatty substances inside the lumen or because of formation of plaques as a result of hemorrhage. As the disease gets progressed, the arteries/arteriole gets constricted. The flow behavior in the stenosed artery is quite different than one in the normal arteries. Also, stresses and resistance to flow are much higher in stenosed arteries in comparison to the normal ones. Having knowledge on flow parameters, such as velocity, flow rate, pressure drop will aid bio-medical engineers in developing bio-medical instruments for treatment (surgical) modalities. Hence fluid mechanics aspect of arterial stenosis have received a considerable attention in the recent past. In view of its importance, a fluid mechanics aspect of arterial stenosis has been undertaken in the present studies. While modeling blood flow in a stenosed tube, it was initially assumed that, the flow obeys Newtonian hypothesis and the flow variables have been computed by using basic NavierStoke's equations [15]. Later, the model has been extended by assuming that, it obeys non-Newtonian hypothesis and developed the model for either Casson fluid or for Power-law fluids [6,7], and showed that, under low shear rates, the model could be best described by this representation. During constriction (stenosis), the lumen of an artery gets considerably reduced thereby, size effects (particle size (mainly red blood cells) to tube diameter ratio) influences the flow characteristics significantly [8,9]. Eringen [10], Cowin [11], and Stokes [12], have proposed micro-continuum theories for accounting size effects in the fluid flows. Later, Ariman et al. [13], Valanis and Sun [14], Chaturani and Pralhad [15] have applied micro-continuum theories for studying the blood flow models. The advantage of application of micro-continuum theories to blood flow model is that, anomalies of blood flows such as Fahraeus and Lindquist effect (FLE) that is dependence of blood viscosity on tube radius, inverse FLE (increase of blood viscosity with decrease in tube radius (15 lm)), and existence of peripheral plasma layer near the tube wall [8,9], can be accounted. In the present model we have accounted micro-continuum fluids proposed by Stokes [12]. These [12] fluids have also been referred as couple stress fluids wherein the parameters , accounts for the size effects in the flow field. Higher value of implies that the flow is tending a g a towards Newtonian whereas, the lower values of implies that the flow has dominance of particle a size effects. is the parameter which accounts for the effect of local viscosity due to particles in g addition to bulk viscosity of the fluid (l). The application of micro-continuum fluids for the stenosed model has been successfully studied by Pralhad and Schultz [16], Sinha and Singh [17] and Srivatsava [18]. The model discussed so far on micro-continuum fluids on stenosis accounts for only viscous effects. Whereas in the flow both the inertia and viscous effects plays an important
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
205
role. In view of the same, an effort has been made in the present model, to account for both inertia and viscous terms by assuming blood as a couple stress fluid.
2. Analysis It is assumed that blood flow is represented by a homogeneous and incompressible couple stress fluid of constant viscosity l, and density q. The constitutive equations and equation of motion for couple stress fluid [12] are Tij;j q dVi dt 1 2 3 4
A eijk Tjk Mji;j 0
Ii;j pdij 2ldij lij 4gxj;i 4g0 xi;j ;
A where Vi the velocity vector, Iij and Tij are the symmetric and antisymmetric part of the stress tensor Tij respectively. Mij is the couple stress tensor, lij is the deviatoric part of Mij , xij is the vorticity vector, dij is the symmetric part of the velocity gradient, g and g0 are the constants associated with the couple stress, p is the pressure and other terms have their normal meaning of the tensor analysis.The appropriate equation from (1)(4) for the two-dimensional flow subject to the additional conditions
Re0
2d ( 1; L0
2R0 % o1 L0
5
together with the equation of continuity can be written as ! ou ou op q v u r2 lu gr2 u or oz oz ! ov ov op q v u r2 lv gr2 v or oz or ou ov v 0; oz or r where 1 o o 2 r r : r dr or
6 7 8
9
Re0 is the Reynolds number, R0 is the tube radius of the unobstructed part (see Fig. 1), ds is the maximum height of the stenosis, L0 is the total length of the stenosis, v and u are the radial and axial velocity components in the r and z directions respectively. The geometry of the stenosis which is assumed to be symmetric is given by
206
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
Fig. 1. Flow geometry of the stenosed tube.
& Rz where
1 s =21 cos2p=L0 fz d L0 =2g; d 1; s ds =R0 : d
d 6 z 6 d L0 ; otherwise
10
Rz Rz=R0 ;
11
Following the order of magnitude analysis of Forrester and Young [19], Eqs. (6) and (7) can be approximated as u ou ou 1 op 1 2 v rw oz or q oz q 12 13
op 0; or where w lu gr2 u: Eq. (12) can now be integrated across the tube to obtain Z Rz o 1 op R2 R ow 2 ru dr ; oz 0 q oz 2 q or R
14
15
where the boundary condition u v 0 at r Rz has been applied. The integrated form of the continuity equation is Z Rz 2pru dr; 16 Q pR2 U
0
where U is the mean velocity at any given cross section with radius Rz, and Q is the volumetric flow rate. We now assume that the radial dependence of the axial velocity can be expressed as a fourth order u 17 Ag1 Bg2 Cg3 Dg4 E 1 1 1 U polynomial of the form, where g1 1 r=R, U is the center line velocity, A through E are, as yet, undetermined coefficients. These coefficients are evaluated from the following boundary conditions:
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
207
i at r R; u 0; where C 1 Ha2 =4I0 a
ii at r 0; ou=or 0;
iii at r 0; u U 18
iv at r R; op=oz r2 w;
v at r 0; o2 u=or2 2U =R2 C ;
g H 21 =a2 =1 f1 I1 a=aI0 ag g pffiffiffiffiffiffiffiffi a a g a Rz ; R0 =l; l g=l; g0 =g:
19
The last initial condition is obtained by the assumption that, at r 0, the velocity profile is not parabolic but of couple stress fluid one. That is u U 1 r=R2 H f1 I0 r=a=I0 ag: 20
So that the second derivative of u with respect to r, at r 0 can be approximated by the condition (v). Here and are couple stress parameters. Eq. (17) reduces to a g u 21 Ag1 6 C 3Ag2 8 2C 3Ag3 3 c Ag4 ; 1 1 1 U where A S1 =S2 a2 k=S2 22 S1 12 2C a2 156 46C S2 7a2 55 dp : dz Eq. (21) can be substituted into Eq. (16) to give ! 30S2 dp 4 2 U 2 Q pR a =15lS2 ; dz pR S3 k R2 =lU where S3 108a2 11a2 C 147C 972: 23 24 25
26
27
With the substitution of Eq. (20) into the left-hand side of Eq. (15), the integral reduces to d 1 dp R2 R ow 2 2 4R U T1 ; 28 dz q dz 2 q or R where 1 H 4H 2H 2 H 2 2 8H D1 D1 D T1 H 2 2 6 2 a a 2 1 a3 and D1 I1 a : I0 a 30 29
208
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
Expressions for U and ouR from Eqs. (16) and (21) are now substituted into Eq. (28) to give an or equation, which when combined with Eq. (26) yields the pressure gradient dp 8qQ2 B1 B2 dR 60lQB3 ; 31 dz dz p2 R5 pR4 where B1 S2 S3 =S5 , B2 Jh 1=3, B3 S2 S4 =a2 S5 . Jh H 2 D2 HB4 D3 B4 =2a3 D4 H=a3 D5 D2 4 aD3 4D2 8=a aD1 1 1 D3 2 4=aD1 D2 1 D4 8a a3 16D1 D5 a3 24a 48D1 8aD2 1 S4 12a 2C a 12a 6C a 288C 288 S5 588a4 63a4 C 810a2 C 5880a2 2475C 9900 B4 21 aa2 1 aD1 2 D6 g g D6 2a 1 fa1 D2 2D1 g: g 1 Eqs. (26) and (31) are now being substituted into Eq. (21) to give the velocity u as a function of r and z. u dR T2 g1 T3 g2 T4 g3 T5 g4 Rz2 T6 g1 T7 g2 T8 g3 T9 g4 Rz3 Re0 B2 ; 1 1 1 1 1 1 dz U0 33 where T2 E1 E2 E1 E3 E4 ; T3 E2 E5 E3 E5 3E4 ; T4 E6 E2 E6 E3 3E4 ; T5 E7 E2 E7 E3 E4 ; T6 E1 E8 E9 ; T7 E5 E8 3E9 ; T8 E6 E8 3E9 ; T9 E7 E8 E9 ; E1 S1 =S2 ; E2 30S2 =S3 ; E3 120S2 S4 =S3 S5 ; E4 a2 =S2 60S2 S4 =a2 S5 ; E7 3 C S1 =S2 ; E5 6 C 3S1 =S2 ; E8 8a2 S2 =S5 ; E6 2C 8 3S1 =S2 ; 34 E9 4a2 S3 =S5 :
4 4 2 2
32
3. Shear stress Shear stress has been computed by two different methods. 3.1. Method-I By using the relation ss Rz dp : 2 dz 35
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
209
Substituting for dp from Eq. (31), and evaluating the shear stress at the maximum height of the dz stenosis (i.e., at z d L0 =2), we get 30lQB3 : 36 pR3 Let s0 represents the shear stress for the Newtonian case (i.e., when a approaches infinity), and in the absence of stenosis (Rz R0 ) Eq. (36) reduces to smax s0 4l1 Q ; pR3 0 15B3 l 3 Rz ; 2 l l=l1 : 37
where l1 is the viscosity of plasma (Newtonian fluid). Using s0 for non-dimensionalizing, we get sR where sR smax =s0 ; 3.2. Method-II By using relation (35), and substituting for dp from Eq. (31), we get dz ! 4 3 sw 4S2 S3 R0 1 dR 60S2 S4 R0 1 Jh 2 2 3 dz S5 Rz a S5 Rz Re0 qU
0
38
39
40
for separation re-attachment point sw2 0, which implies qU 0 dR R0 15S4 ; Re0 2 1 dz R a S3 3 Jh where Rz & R0 R0 ; ds
2
40a
1 cospz=z0 ;
2z0 L0 : otherwise
41
4. Resistance to flow Resistance to flow, which is of physiological importance, has also been computed by two methods. 4.1. Method-I The resistance to flow, k is defined by P0 P k ; Q
42
210
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
where P0 and P is the pressure at the entry and exit level respectively. Integrating Eq. (31) using the condition that p p0 at z 0 and p P at z L, we have Z L Z 8qQ2 L B1 B2 dR 60lQ B3 P0 P 2 dz: 43 dz 5 dz p p R R4 0 0 Using Eq. (10) for Rz, for the substitution of dR in Eq. (43), k is given by dz Z dL0 60lL 60lB3 8qQB1 B2 dR 1 L0 B3 dz; k pR4 pR4 p2 R5 dz d 0
44
where B B3 at Rz R0 , and L0 L0 =L. 3 Let k0 be the resistance to flow for the Newtonian case in absence of stenosis. Eq. (44) reduces to k0 8l1 L : pR4 0 45
Using k0 for non-dimensionalization with k, Eq. (44) simplifies to Z p Z d 15B l 3 15L0 B3 Re0 ls R0 p B1 B2 sin h 1 L0 l kR dh dh; 4 2 4p 4 Rz5 p Rz p where k kR ; k0 4.2. Method-II Using Rz in the form of Eq. (41) and simplifying the integration term, we get P0 P
2 qU 0
46
2R0 qU 0 Re0 ; l
2p zd h L0
L0 2
! ;
R0
R0 : L
47
2B1 B2 Rz4 1
40B3 Z 0 A; A1 p Re0
48
where A
3 5 2 2:5 2 f3 2 0:5 g sin h dA2 dA2 dA2 d 4A2 2A2 Rz 2Rz 4A1 Rz 1 1 2 pffiffiffiffiffi 6A3 2:25 A2 h d pi 1 2 pffiffiffiffiffi tan ftanh=2= A1 g 2 A5 1
sin h d
3
5A2 sin h d
2
49
d; A1 1
A2 d=2;
1 ds ; d R0
R0 l
Z0
Z0 ; R0
U0
Q : pR2 0
Here `a' is assumed to be of the form models.
in order to compare with the other theoretical
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
211
5. Comparison to other theoretical models The present model has been compared with the other theoretical models [17,19]. The present model can be reduced to the model of Forrester and Young [19] by allowing the parameter a a Rz tending to infinity (which implies the absence of size effects) in the equation for velocity (Eq. (33)), shear stress (Eq. (40)), and in the equation for the resistance to flow (Eq. (48)). The reduced form is given below. ! 2 u 308 1204 2 4 3 4 4 dR 2 22g1 g1 1=Rz 50 g g g g Re0 1=Rz3 : 1575 1 1575 1 5 1 15 1 dz U0 sw 616 8 4 dR 1=Rz3 : 51 1=Rz 2 1575 dz Re0 qU
0
P0 P qU 0
2
4 194 z 1=Rz 1 160 =3p Re0 A1 A : 225
52
The results of the present investigation have also been compared with the analysis of Sinha and Singh [17]. The results of their investigations have been modified as per the reasons stated in Ref. [20] and the modified results are given below. s R k R l Rz3 B5 : L0 l 2p Z
p
53 dh Rz4 B5
l1
3L0 =2B 5
: s s=s0 ; R k k=k0 : R
54 55
p
B5 1 2H=aa 2D1 ;
B B5 jRzR0 ; 5
While comparing the results of present model with that of Forrester and Young [19], the equation of the tube radius of the present model (Eq. (10)) and that of Forrester and Young [19] (Eq. (41)) have been used independently and while plotting the graph (Figs. 68) these changes have been accounted. The h variable which is mentioned in Eq. (47) is a part of the substitution for the computation of the integral in Eq. (46) by the Simpson's rule. It should not be treated as another form of Rz rather it is the same form as referred in Eq. (10).
Table 1 Experimental values of rotational viscosity, lR , for different haematocrit [13] Haematocrit (%) 10 20 40 lR (cP) 0.78 0.92 0.98
212
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
Table 2 Experimental values of l (from [21]), and theoretical values of , and g (from [14]) for different diseases a g Diseases Normal blood Polycythemia Plasma cell dyscrasias Hb SS (sickle cell) Viscosity (l, cP) 3.81 6.75 4.99 3.29 a 2.63 3.50 3.01 2.44 g )0.5 )0.5 )0.5 )0.5 g (cP cm2 ) 0.5 0.5 0.5 0.5
Fig. 2. Variation of shear stress (sR ) with stenosis height (s ) for different , and concentration. d a g
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
213
6. Computation details The flow variables have been computed for the different suspension concentration of haematocrits (10%, 20%, and for 40%) and for the diseases such as, polycythemia, sickle cell (Hb SS), and for plasma cell dyscrasias. The results have also been compared to the case of normal blood. The value of l, when computing for different suspension concentrations, is chosen to be of the form l l1 lR where l1 is the viscosity of plasma (Newtonian fluid) and lR is the rotational viscosity due to suspension of blood cells. The value of l1 is assumed to 1.2 cP (centi Poise) and lR for different suspension concentrations are taken from Ariman et al. [13] and are shown in Table 1. For the computation of axial velocity, resistance to flow and the shear stress distribution for the different diseases, the values of l have been taken from Chien [21]. The value of has been fixed to g )0.5 for the computation of flow parameters for all the diseases and for the normal blood (Table 2). Resistance to flow for different suspension concentrations (Eqs. (46) and (54)) have been computed by using Simpson's rule. The values of s and L0 are varied from 0 to 0.5, and 0.5 to 1.0 d respectively. Size effects parameters and are varied from 2 to 20 and )1.0 to 1.0 respectively. a g Lower the value of more is size effects pronounced in the model, whereas the higher values of a a accounts that the model is tending towards the Newtonian results. values are restricted between g )1 and 1 only (Ref. [14], Eq. (41)). Chaturani and Upadhya [22], and Chaturani and Pralhad [23] have discussed in detail on the aspect of choice of parametric values of , for the blood flow a g
Fig. 3. Variation of shear (sR ) stress with for different concentrations (s 0:20). a g d
214
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
modeling. The values of R0 and 0 , have been taken to be 0.02 and 4 respectively. The values of z a z=z0 have been varied from )1 to +1. It is worthwhile to mention at this stage that, though z and are related through l, however their values are chosen separately due to the non-availability g of the experimental values on g and g0 . The reasons for computing flow variables for the different suspension concentrations and for different diseases are due to the fact that the model can be used for general-purpose reference. Once the volume concentration of blood cells (haematocrit) for a particular disease is known, the flow variables can be computed by making use of Tables 1 and 2 and the method of interpolation and extrapolation. In comparison of the flow variables, the tube geometry is assumed to be of the size 0.95 cm. The reason for accounting higher tube size is for comparison with the existing Newtonian models [19]. In doing so, it is assumed that the flow is steady and laminar. The assumption is quite acceptable for the reason that, the Reynolds numbers chosen for the computation purpose is ranging from 10
Fig. 4. Variation of resistance to flow (kR ) with stenosis height (ds ) for different , L0 and concentrations ( 3:0, g a Re0 10).
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
215
to 400 only. Once the results are confirmed (after comparing with the other theoretical models) that the present model yields satisfactory results, the model then can be computed for the lesser tube diameter also.
7. Results The results on shear stress (Eqs. (38) and (53)) and on resistance to flow (Eqs. (46) and (54)) for different suspension have been shown in Figs. 25. It is observed that the values of shear stress increases with the increase of stenosis height and decrease with the increase of couple stress parameters and . Also, it indicates that the values of the present model are lower in comparison a g to the values of Sinha and Singh [17] for the parametric values of > 3 but for 6 2, the values of a a the present model are higher in comparison to Ref. [17]. The results on axial velocity (Eqs. (33) and (50)), shear stress (Eqs. (40) and (51)) and resistance to flow (Eqs. (48) and (52)) for different diseases are shown in Figs. 68. The results on shear stress
Fig. 5. Variation of resistance to flow (kR ) with for different and concentrations (L0 0:5, and s 0:20). a g d
216
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
u r Fig. 6. Variation of axial velocity U with tube radius R at different axial locations (z), for different Reynolds 0 numbers (Re0 ) and diseases (() recent model; (- - -) Forrester and Young [19]).
and resistance to flow agree well with the quantitative observation of Forrester and Young [19] (shown as dotted lines). The values are lower in case of shear stress distribution, and are higher in case of resistance to flow in comparison to Forrester and Young [19]. The comparison between the present model and that of Sinha and Singh [17] is rather poor (except for the trends) for the set of values computed for the shear stress and the resistance to flow as shown in Figs. 25. Whereas the comparison between the present model and that of Forrester and Young [19] are quite good (Figs. 68). This is for the reason that the present model, and that of Forrester and Young [19], accounts
R.N. Pralhad, D.H. Schultz / Mathematical Biosciences 190 (2004) 203220
217
Fig. 7. Variation of Shear stress sR2 with axial location Zz0 for different Reynolds numbers (Re0 ) and diseases (()
qU 0
present model; (- - -) Forrester and Young [19]).
for both the inertia and the viscous terms whereas Sinha and Singh [17] model accounts for only viscous terms. The values of axial velocities appear to be higher in comparison to Ref. [19] for lower values of Reynolds number 6 100, but for higher values of Reynolds number, Newtonian values predominate over the present ...