0021-925ap1
J. Ilimdda Vol. 24, No. 11, pp. 106%1077. 1991. Rirdd in Chat Briti
MATHEMATICAL MODELLING OF FLOW THROUGH IRREGULAR ARTERIAL STENOSIS PETER R.JOHNSTON
s3.al+.w
Pelgamoo Flu8
plc
AN
and DAVID KILPATRICK
Department of Medicine, University of Tasmania, 43 Collins St, Hobart 7000, Australia Abstract-A mathematical model of flow through an irregular arterial stenosis is developed. The model is two-dimensional and axi-symmetric with the stenosis outline obtained from a three-dimensional casting of a mildly st&osed artery. Agreement between modelled and experimental pressure drops (obtained from an axi-symm&ic machined stenosis with the same proBe) is excellent. Resultsiare also obtained for a smooth stenosis model, similar to that used for most mathematical modelling studies. This model overestimates the pressure drop across the stenosis, as weIf as the wag shear stress and~separationReynolds number. Also, the smooth model predicts one instead of three recirculation zones present in the irreguhtrmodel. stenosis is modified to increase the severity from 48 and 87% areal occlusion, while g the same general shape. This has the effect of increasiug the pressure drop by an order of decreasing the number of recirculation zones to one, with a lower separation Reynolds number.
1. INTRODUCTION
Occlusive arterial disease is the major cause of death in western society. The common form of arterial narrowing (stenosis) is that caused by atheroma, a deposition of fats and fibrous tissue in the arterial wall. The effect of this is to reduce arterial lumen to an extent which reduces blood delivery to the region beyond the narrowing. Several modelling studies have been performed on a single stenosis. The solution methods for these have varied from the momentum integral method (obtained by neglecting the radial pressure drop) (Forrester and Young, 1970; Mbrgan and Young, 1974) to sophisticated numerical methods (Lee and Fung, 1970; Lkshpande et al., 1976; Johnston and Kilpatrick, 1990). These studies have all assumed that the stenosis could be represented by a smooth mathematical function, for example, the cosine curve. In reality, this is not generally the case. Arterial stenoses, whilst having the general trend of a smooth curve, contain many small valleys and ridges, similar to a mountain range. To study this ,problem in greater detail, published data (Back et al., 1984) were used to define the outline of the stenosis. This was collected from a casting of a I& circumflex coronary artery with mild, ditTuse atherosclerotic disease. The original purpose of the casting was to observe the flow from an experimental point of view (Back et al, 1984). Here it is used to develop an axi-symmetric mathematical model for the flow. The model is validated by comparing pressure drop predictions with those measured in experiments using a machin axi-symmetric stenosis. With the vali$3 ted model, a detailed study of the flow patterns can be undertaken. Velocity and wall shear stresses are considered and these show that Receiwd infinaljii
6 May 1991.
several regions of recirculating flow are present at high dow rates. These flow patterns are compared with a smooth stenosis of the same length and severity. Also, the stenosis is made more severe by decreasing the quoted radius values. In addition, the pressure drops, velocity fields and wall shear stresses are studied and in turn compared with a smooth stenosis, as indicated above. 2. STENOSISMODEL As mentioned above, the profile of a left circumflex coronary artery was obtained in the following manner (Back et al., 1984). Liquid silicone rubber was injected into the artery at physiological pressure and a polyurethane negative was taken from the hardened sihcone rubber casting. The cross-sectional area of the negative cast was computed at 50 micron intervals from two orthogonal X-ray images taken along the cast axis. A solid axi-symmetric model, with radii calculated on the assumption of circular cross-sea tional areas of the original section, was fabricated by machining an ahtminium rod to the appropriate diameter and then making a polyurethane casting. The numer&al model of the stenosis was constructed from these data Table 1 shows the variation in arterial radius along the axis of the stenosis. Here, the stenosis severity is 48% area1 occlusion. The stenosis is also shown diagramaticahy in Fig. 1. The original use for the casting was to measure the pressure drop across a three-dimensional stenosed arterial section (Back et al., 1984). Further pressure drop experiments were performed with the axi-symmetric casting. It is these results which are used for comparison with and validation of the model in this paper. Other physical measurements required were also taken from the literature (Back et ai., 1984). The overall length of the stenosis, including short inlet and
1069
1070
P. R. JOHNSTONand D. KILPATIUCK Table 1. Variation of stenosis radius along the axis Axial Position 0.100 Z:E 0.942 1.039 1.136 1.201 1.299 1.396 1.494 1.591 1.688 1.786 1.883 1.981 2078 2.143 2.240 2.338 2.435 2.532 2.630 2.727 2.825 2922 3.019 3.084 3.182 3.279 3.377 3.474 3.571 3.669
Axial Position
Radius
Axial Position
Radius
3.766 3.864 3.961 4.058 4.123 4.22 1 4.318 4.416 4.513 4.610 4.708 4.805 4.903 5s00 5.065 5.162 5.260 5.357 5.455 5.552 5.649 5.747 5.844 5.942 6.006 6.104 6.201 6.299 6.396 6.494 6.591 6.688 6.786
0.426 0.422 0.418 0.414
6.851
0.410 0.422 0.433 0.441 0.452 0.449 0.452 0.452 0.449 0.441 0.433 0.433 0.441 0.449 0.463 0.470 0.474 0.478 0.474 0.478 0.48 1 4.492 0.492 0.488 0.485 0.481 0.485 0.495 0.492 0.488 0.500 0.500 0.500
0.500 0.500 0.500 0.498 0.502 0.495 0.505 0.502 0.502 0.502 0.505 0.508 0.508 0.505 0.505 0.502 0.505 0.505 0.498 0.495 0.488 0.488 0.485 0.474 0.478 0.470 0.452 0.452 0.441 0.441 0.437 0.437 0.433
2 4 Dimen8ionlear Pig. 1. Pro!&
of irrc~r
outlet flow tubes, was 3.246 cm and the internal diameter of the inlet flow tube (Do) was 0.308 cm. This diameter was used to render the model non-dimensional, giving the Reynolds number as Re=-
U,D,P P
EE 7.143 7.240 7.338 7.435 7.532 7.630 7.727 7.825 7.922 7.987 8.084 8.182 8.279 8.377 8.474 8.571 8.669 8.766 8.864 8.929 9.026 9.123 9.221 9.318 9.416 9.513 9.610 9.740 10.000 10.300
:z Ok6 0.401 0.401 0.406 0.410 0.414 0.410 0.406 0.401 0.397 0.389 0.384 0.371 0.362 0.366 0.366 0.366 0.362 0.362 0.366 0.366 0.371 0.371 0.371 0.371 0.384 0.401
6 Axial
8
10
Position
stenosis (Back et I& 1984). and the dimensionless pressure as PC-L pu:/2 where u,, is the average axial velocity in the unoccluded tube, p the fluid density, p the viscosity and p the pressure.
Mathematical modelling of an irregular arterial stenosis 3. GOVERNINGEQUATIONS
Blood flow through an axi-symmetric stenosis can be modelled in two dimensions using a cylindrical coordinate system. The dimensionless governing continuity and Navier-Stokes equations for this system are:
au, V, au,
a,+-+-“o, r
au, V, -jy
au, ap Y&=
--g
I +-
au, v, F+“’
az
au,
a3 I au --L+_-.L_‘+:
[ a9 ap I
Re
a5
ar r2
I
a%,
-=-jg+z a2
v
I
a221 au,
-g+;dr+j$[
,
(2)
a%, 1
’
1071
used (191 in the longitudinal direction and 19 radially), with a high concentration of elements near the wall and in the region of the stenosis (approximately 0.3 by 0.2 mm), decreasing with distance away from the stenosis. A system of algebraic equations is obtained from the mesh points and element connectivity by integrating the shape functions from the Galerkin weighted residual representation of the governing equations. The pressure degree of freedom has been eliminated using a penalty function formulation (parameter, E- 1 x 10T6) and is recovered after the velocities have been calculated. This has the effect of reducing the size of the system of equations. The system of equations is solved iteratively using up to three successive substitutions followed by quasi-Newton iterates with Boyden’s update. Generally, only five or six iterations were required for the solution to converge, leaving residuals of about 10e4 over the whole grid.
(3) where r and z are the dimensionless coordinates, scaled with respect to D,, with the z-axis located along the symmetry axis of the artery. No secondary or swirling flows have been allowed so the total velocity is defined by the dimensionless radial and axial components, u, and v,, scaled with respect to u,. These equations assume that the fluid is incompressible and Newtonian, which makes it reasonable to compare with the experimental results as water was the fluid used. Unfortunately, this theory does not lend itself to accurately modelling blood flow. The velocity boundary conditions on the arterial wall are the usual no slip conditions at r=R(z).
v,=vz=O
(4)
Flow upstream Ifrom the stenosis was assumed to have a velocity profile corresponding to Hagan-Poiseuille flow througb a long circular tube of constant crosssection v,=2
r2
( ) l--
,
R,2
v,=o,
(6)
(R, = D,,/Z). It was decided, after numerical study, to apply this boundary condBion eight arterial radii upstream from the beginning of the stenosis, thus avoiding interaction between the stenosis and the entering parabolic velocity profile. Downstream, the velocity is left free but conditions of &ro normal and tangential stress are applied. These ~conditions were applied eighty radii downstream from the beginning of the stenosis in order that fully developed flow be reestablished (Engtsman, 19% The govern+8 equations, being highly nonlinear, were solved nu+neri&ly using the FIDAP computational 5uid d ‘amice package, utilising the Gale&in bod (En&man, 1990,1982). A non5nite element F uniform grid of 13429linear quadrilateral elements was B?I 24:11-G
4. MODEL VALIDATION The model was validated by comparing modelled pressure drops with the experimental values given in the literature. Pressures are calculated by differencing the axial pressures at z=O.42 (0.13 cm) and z- 10.118 (3.12 cm) from the beginning of the stenosis. Note that these points are actually in the solid tubes attached at each end of the irregular section of the stenosis. A plot of the comparison is shown in Fig. 2. Good agreement is obtained for most of the range of Reynolds numbers considered. The difference is only of the order of a few per cent. Larger errors are observed with increasing Re and are due probably to the occurrence of 5ow instabilities or possible transition to turbulence (Deshpande et al., 1976), neither of which is accounted for in the model. These are regions where the 5ow patterns deviate from the smooth laminae assumed in the model. Over the physiological flow range (Re = 100400) the agreement is very good. This good agreement is the basis for using the model for a detailed study of the 5ow patterns obtained in the stenosis. Despite the fact that differences between the experiment and the model occur at higher Reynolds numbers, results are presented at these higher values for comparison with previous studies (Deshpande et al., 1976; For-rester and Young, 1970; Morgan and Young 1974).
5. SEPARATlONAND RFXIRCULATINCFLOW Points of flow separation are observed when the normal wall shear stress changes from positive to negative. (DetaiR of normal wall shear stress comput&ions arc given in Engleman, 1990.) The flow is said to reattach when the stress becomes positive again. A normalii wall shear stress is obtained by dividing calculated values by the value in the un-
P.R. iOHNSTON and D. KILPATIUCK
1072
II)
100 Reynolds Number
J 1000
Fig.2. Comparisonofmodelled and experimental pressure drops.
I
8,
I 11
7 -
I’
I' 11 ;I
6 -
Re-20 Re=lOOO ---- _
I’ I 1III !/I1
5 4 -
n
-1
’
0
I
2 4 Dimensionless
6 8 Axial Position
10
Fig.3.Wallshearstrewprolila along the stenosis at Reynolds numbers of 20 and 1000.
occluded tube section, far upstream from the stenosis (i.e. the value obtained from Poiseuillean flow). In the stenosis under consideration, there is no flow separation observed at low Reynolds numbers. At a Reynolds number of 600 a small region of recirculating flow develops between Z= 7.12 (2.19 cm) and z = 7.24 (2.23 cm) downstream from the beginning of the stenosis. Thii represents the separation Reynolds number, that is, the Reynolds number at which flow separation is first observed. When the Reynolds number is increased to 800, two recirculation zones are observed. The 6rst is a larger region in the same position as at Re=600, and the second is further downstream between z= 8.82 (272 cm) and z= 8.88 (2.74 cm). Increasing the Reynolds number further to 1000 produces another recirculation zone which is between the two at Re = 800, from z = 8.21(2.53 cm) to
z = 8.54 (2.63 cm) downstream from the beginning of the stenosis. It is interesting to locate the separation and reattachment points on the stenosis curve (refer to Fig. 1). The first region of flow reversal is in the large valley after the narrowest point of the stcnosia The second and third regions are in the next two valleys beyond the first. However, there is no recirculation in the fourth valley, near the exit of the stenosis. It is possible that a region may develop in this area at a higher Reynolds number. Figure 3 shows the wall shear stress vs axial position at Reynolds numbers of 20 and 1ooO. Note that the two curves are very similar in shape. Although the figure shows plots for only the extreme values of the Reynolds numbers considered, all curves for the intermediate Reynolds numbers lie between these two. The
Mathematical modelling of an irregular arterial stenosis
1073
where z,, is the half-length, 6 the width and z1 the centre of the stenosis. Here thii is used for comparison with results obtained from the irreguhu stenosis The centre of the smooth stenosis is taken at the narrowest point of the irregular one. Figure 4 shows a comparison of the smooth and irregular s&noses. A mesh containing 3629 linear quadrilateral elements, of similar size to those utiri above, was used to solve this model. 6. COMPARISONWllIi A SMOOTHSl’ENOSIg Comparing pressure drops over the two stenoses shows that the smooth stenosis predicts higher values Most model&g studies performed use a smooth than the irregular stenosis (Fig 5). The two curves are stenosis of the form parallel but the model of the irregular steno& gives much better agreement with the experimental results. A plot of the wall shear stress shows that the smooth (z)= stenosis predicts a much narrower range of values (Fig. 6). In the converging section of the stenosis the z>hJ, (7) irregular model predicts lower shear stresses, but at
points of separation and reattachment are clearly visible as the crossing points of the horizontal axis. Maximum shear stress occurs at the narrowest point in the stenosis. The expansion downstream of the narrowest point results in a rapid drop in shear stress in that region. In fact, the shear stress profile retlects very closely the outline of the stenosis.
l+u)s(y)), i:: loo0 (compared with 600 for the irregular model). It is also further away from the narrow point of stenosis than seen with the irregular stenosis. Another point of comparison is the axial velocity. Figure 7 shows that velocities through the irregular stenosis are less than those in the smooth one. However, it appears that the shapes of the velocity profiles are similar, implying that the tendency for flattened velocity proflles through the stenosis is the same in both models.
2
,4
7. MORE SEVEW SIXNOSES
To study the effects of a more severe stenosis, the original stenosis profile was modified to achieve a severity of 87% areal occlusion, while maintaining the same general shape. The results profile is shown in Fig 8, along with the original profile. A plot of wall shear stress is shown in Fig 9 for Reynolds numbers of 20 and 1000. These curves are smoother than those for the l& severe stenosis (Fig 3), indicating that the more severe stenosis tends to appear smoother in the diverging section. Also, the values are an order of magnitude larger than for the lesser constriction. There is only one region of flow recirculation which encompasses the three regions observed previously and it occurs for Reynolds numbers > 60. The presence of only one region is due to the
-8
6
Dimensionless Axial Position Fig. 6. Wall shear stress pro&s for irrqular and smooth models at Reynolds numbemof 20 and 1ooO.
‘\ ‘\ \
‘\ \
, I L
@
\\
\
\\
\
\
\
\
\
\
\
.35 .l .15 .2 .25 .3 Dimensionless Radial Position
0.05
.4
Mathematical modelling of an irregular arterial stenosis
1075
.s .4
.3
.2
Original Modified
.l
0
2
4
Dimensionless
Axial
6
Data Data
---_
8
10
Position
Fig. 8. Modified irregular stenosisand the original stenosis.
80
Re=lOOO ----
70 60 50 40 30 20 10 0 A”
0
2
4
Dimensionless
8
6 Axial
10
Position
Fig. 9. Wall shear stress profiles along the stenosis for the mcxlilkd irregular stenosis at Reynolds numbers of 20 and 1000.
Modified Smooth 0
’t 0
Original Data Approximation
----
1 I
2
4
Dimensionless
6 Axial
8 Position
Fig. 10. Smooth approximation to the mod&d original data
10
P.R.JOHNSON andD. KILPATRICK
1076
higher velocities achieved in the stenosis, allowing little opportunity for the flow to reattach to the walls even when there is a slight recontraction in the pro&. As previously stated, a smooth stenosis profile was developed to compare with this more severe stenosis. Figure 10 shows the two profiles. Rressure drop comparisons are shown in Fig 11. The two profiles are parallel again, as for the less severe stenosis, but the smooth profile model again overestimates the total pressure drop, this time to a greater extent. A plot of wall shear stress at Reynolds numbers of 20 and loo0 is shown in Fig. 12. Again the smooth profile model overestimates the stress in the converging section but underestimates it at the narrowest point. In the diverging section of the stenosis, the smooth model matches the irregular model very closely and, more importantly, gives a very good prediction of the point of flow separation. The separa-
tion Reynolds number for this model is < 60, which is an underestimate of the irregular model value. This is in contrast to the less severe situation, where the smooth model overestimated the separation Reynolds number. These observations are due probably to the similarity in stenosis pro&s in the diverging section. Finally, Fig. 13 shows axial velocity pro&s for the two severe stenoses. The smooth profile predicts higher velocities, as was observed in the less severe case. However, there is a difierence in the shape of the velocity profiles. The typical blunt velocity profile of flow through a stenosis is more exaggerated with the irregular profile. a
CONCLUSION
The mathematical model of flow through an irregular stenosis presented here shows excellent agreement
1000
Modified Irregular Model I Smooth Stenosis ---- ’
1.00:
.A”
~
10
100 Reynolds Number
Fig 11.Prcssuredropformodificdirrcgularsnd ~00th mod&.
no-20
g
g
ao-
lb-20
70.
ti
t 60?z 50 2 g 40. 3
(Irregulars
It,-1000 (Irregulars
It.4000
,swoth) @mothl
---1 I, ....'...I,'I ----
II
30.
fE 20%A 01 10 4
ca
-10 ' 0
2. 4 6 8 Dimensionless Axial Position
Fig. 12. Wall shear strcsa profrlcsfor modi!icd irregular nod smaoth mod& at Rcynold8 numbers of20 and1000.
Mathematical modelling of an irregular arterial stenosis
------------t
--__
1077
Irregular Model *_._\ Smooth Model ----
6
0 0.020.040.060.08 .l .12 .14 .16 .18 .2 Dimensionless Radial Position Fig. 13. Axial velocity profiles for the modified irregular and smooth models at z=6 and Re=200
with publish+ experimental data for the same stenosis. This made possible a detailed study of flow patterns inthkstenosis and the identification of three recirculation .xones in the diverging section of the stenosis. A comparison of results from the irregular model with the more common smooth model showed that the latter overlooked several fluid dynamical aspects of the irregular situation. These were overestimates for pressure droR, wall shear stress (except at the maximum) and separation Reynolds number. Also, only one recirculation zone was identified. These errors are due mainly to the simple shape assumed for the stenosis. When a more severe stenosis of similar shape was studied, only one recirculation zone was observed, due to the higher fluid velocity in the narrow section of the stenosis. A smooth stenosis of the same severity again overestimated the pressure drop and wall shear stress but underestimated the separation Reynolds number. It can be concluded from this study that it is necessary to have a good representation of the stenosis profile to make accurate predictions of the fluid flow in the stenosis. This is true especially when looking at recirculation xones. These are postulated to be prime areas for further deposition of atherosclerotic plaques which can lead to a worsening of the occlusion (Asakura and Karino, 1990). Methods based on fibre optic laser Doppler anemometry (Webber, 1982) and ultrasound (Kitney et al., 1990) are available to obtain estimates of the internal protIle of a s&nosed artery. These data could be used easily in a model, such as has been developed here, to give accurate predictions of the art&al fluid flow (assuming a gid artery with a Newtonian fluid) and therefore hi 4 light any potential areas of further deposition or extreme wall shear stress. Clearly, a threedimensional Buid flow model would give even better
insight, but this is not practical from a computational point of view at this stage. Other factors which may be important and should be studied are the inclusion of arterial wall motion and a non-Newtonian model for blood. Acknowledgement-PRJ was ,supported by the National Heart Foundation of Australia. REFRRRNCES
Asakura, T. and Karino, T. (1990) Flow patterns and spatial distribution of atherosclerotic lesions in human coronary arteries. Circulation Res. f&1054-1066. Back, L., Cho, Y, Crawford, D. and Culfel, R. (1984) Effect of mild atherosclerosis on 5ow resistance in a coronary artery casting of man. J. bionw?ck.Engng 1% 48-53. Deshpande, M., Giddens, D. and Mabon, R. (1976) Steady flow through modelkd vascular stenoses. J. Bimechics 9, 165-174. En&man, M. S. (1982) FIDAP-a fluid dynamics analysis &age. Ado. &I@ S&V. 4,163-166. _ _ Enakman M. S. 11990)FIDAP User’s Guide. Fluid Dvnamh International, Evanston, Illinois. Version 5.0 E&t. Forrester, J. H. and Young, D. F. (1970) Flow through a converging-diverging tube and its implications in occlusive vascular dVIZ theoretical development. J. Biomecka?tics3,297-305. Johnston, P. R. and Kilpattick. D. (1990) A mathematical model of paired arterial stenoses. In Computers in Curdiology (Edited by Ripky, K), pp. 229-232. IEEE Computer Society pnsS Chicago. Kit&y, R. I., Straughan, Kz Moura, L., Burr& C. and Rothman M. T. (1990) 3-d visuahsation for the studv of arterial disease and tissue chamcterimtion. In Com&ters in Cardiology(Edited by Ripley, K.), pp. 3-6. IEEE Computer Society Preq Chicago. Lee, J.-S. and Fung, Y.-C. (1970) Flow in looally constricted tubes at low Reynolds numbers. J. appl. MvcA 37,9-16. Morgan, B. and Young, D. (1974) An integral method for the analysis of flow in arterial stenoses. Bull. math. Biol. 36, 39-53. Webber, S. (1982) The haemodynamks of human coronary artery stenosis. M.Sc. thesis, University of Tasmania.
J. Ilimdda Vol. 24, No. 11, pp. 106%1077. 1991. Rirdd in Chat Briti
MATHEMATICAL MODELLING OF FLOW THROUGH IRREGULAR ARTERIAL STENOSIS PETER R.JOHNSTON
s3.al+.w
Pelgamoo Flu8
plc
AN
and DAVID KILPATRICK
Department of Medicine, University of Tasmania, 43 Collins St, Hobart 7000, Australia Abstract-A mathematical model of flow through an irregular arterial stenosis is developed. The model is two-dimensional and axi-symmetric with the stenosis outline obtained from a three-dimensional casting of a mildly st&osed artery. Agreement between modelled and experimental pressure drops (obtained from an axi-symm&ic machined stenosis with the same proBe) is excellent. Resultsiare also obtained for a smooth stenosis model, similar to that used for most mathematical modelling studies. This model overestimates the pressure drop across the stenosis, as weIf as the wag shear stress and~separationReynolds number. Also, the smooth model predicts one instead of three recirculation zones present in the irreguhtrmodel. stenosis is modified to increase the severity from 48 and 87% areal occlusion, while g the same general shape. This has the effect of increasiug the pressure drop by an order of decreasing the number of recirculation zones to one, with a lower separation Reynolds number.
1. INTRODUCTION
Occlusive arterial disease is the major cause of death in western society. The common form of arterial narrowing (stenosis) is that caused by atheroma, a deposition of fats and fibrous tissue in the arterial wall. The effect of this is to reduce arterial lumen to an extent which reduces blood delivery to the region beyond the narrowing. Several modelling studies have been performed on a single stenosis. The solution methods for these have varied from the momentum integral method (obtained by neglecting the radial pressure drop) (Forrester and Young, 1970; Mbrgan and Young, 1974) to sophisticated numerical methods (Lee and Fung, 1970; Lkshpande et al., 1976; Johnston and Kilpatrick, 1990). These studies have all assumed that the stenosis could be represented by a smooth mathematical function, for example, the cosine curve. In reality, this is not generally the case. Arterial stenoses, whilst having the general trend of a smooth curve, contain many small valleys and ridges, similar to a mountain range. To study this ,problem in greater detail, published data (Back et al., 1984) were used to define the outline of the stenosis. This was collected from a casting of a I& circumflex coronary artery with mild, ditTuse atherosclerotic disease. The original purpose of the casting was to observe the flow from an experimental point of view (Back et al, 1984). Here it is used to develop an axi-symmetric mathematical model for the flow. The model is validated by comparing pressure drop predictions with those measured in experiments using a machin axi-symmetric stenosis. With the vali$3 ted model, a detailed study of the flow patterns can be undertaken. Velocity and wall shear stresses are considered and these show that Receiwd infinaljii
6 May 1991.
several regions of recirculating flow are present at high dow rates. These flow patterns are compared with a smooth stenosis of the same length and severity. Also, the stenosis is made more severe by decreasing the quoted radius values. In addition, the pressure drops, velocity fields and wall shear stresses are studied and in turn compared with a smooth stenosis, as indicated above. 2. STENOSISMODEL As mentioned above, the profile of a left circumflex coronary artery was obtained in the following manner (Back et al., 1984). Liquid silicone rubber was injected into the artery at physiological pressure and a polyurethane negative was taken from the hardened sihcone rubber casting. The cross-sectional area of the negative cast was computed at 50 micron intervals from two orthogonal X-ray images taken along the cast axis. A solid axi-symmetric model, with radii calculated on the assumption of circular cross-sea tional areas of the original section, was fabricated by machining an ahtminium rod to the appropriate diameter and then making a polyurethane casting. The numer&al model of the stenosis was constructed from these data Table 1 shows the variation in arterial radius along the axis of the stenosis. Here, the stenosis severity is 48% area1 occlusion. The stenosis is also shown diagramaticahy in Fig. 1. The original use for the casting was to measure the pressure drop across a three-dimensional stenosed arterial section (Back et al., 1984). Further pressure drop experiments were performed with the axi-symmetric casting. It is these results which are used for comparison with and validation of the model in this paper. Other physical measurements required were also taken from the literature (Back et ai., 1984). The overall length of the stenosis, including short inlet and
1069
1070
P. R. JOHNSTONand D. KILPATIUCK Table 1. Variation of stenosis radius along the axis Axial Position 0.100 Z:E 0.942 1.039 1.136 1.201 1.299 1.396 1.494 1.591 1.688 1.786 1.883 1.981 2078 2.143 2.240 2.338 2.435 2.532 2.630 2.727 2.825 2922 3.019 3.084 3.182 3.279 3.377 3.474 3.571 3.669
Axial Position
Radius
Axial Position
Radius
3.766 3.864 3.961 4.058 4.123 4.22 1 4.318 4.416 4.513 4.610 4.708 4.805 4.903 5s00 5.065 5.162 5.260 5.357 5.455 5.552 5.649 5.747 5.844 5.942 6.006 6.104 6.201 6.299 6.396 6.494 6.591 6.688 6.786
0.426 0.422 0.418 0.414
6.851
0.410 0.422 0.433 0.441 0.452 0.449 0.452 0.452 0.449 0.441 0.433 0.433 0.441 0.449 0.463 0.470 0.474 0.478 0.474 0.478 0.48 1 4.492 0.492 0.488 0.485 0.481 0.485 0.495 0.492 0.488 0.500 0.500 0.500
0.500 0.500 0.500 0.498 0.502 0.495 0.505 0.502 0.502 0.502 0.505 0.508 0.508 0.505 0.505 0.502 0.505 0.505 0.498 0.495 0.488 0.488 0.485 0.474 0.478 0.470 0.452 0.452 0.441 0.441 0.437 0.437 0.433
2 4 Dimen8ionlear Pig. 1. Pro!&
of irrc~r
outlet flow tubes, was 3.246 cm and the internal diameter of the inlet flow tube (Do) was 0.308 cm. This diameter was used to render the model non-dimensional, giving the Reynolds number as Re=-
U,D,P P
EE 7.143 7.240 7.338 7.435 7.532 7.630 7.727 7.825 7.922 7.987 8.084 8.182 8.279 8.377 8.474 8.571 8.669 8.766 8.864 8.929 9.026 9.123 9.221 9.318 9.416 9.513 9.610 9.740 10.000 10.300
:z Ok6 0.401 0.401 0.406 0.410 0.414 0.410 0.406 0.401 0.397 0.389 0.384 0.371 0.362 0.366 0.366 0.366 0.362 0.362 0.366 0.366 0.371 0.371 0.371 0.371 0.384 0.401
6 Axial
8
10
Position
stenosis (Back et I& 1984). and the dimensionless pressure as PC-L pu:/2 where u,, is the average axial velocity in the unoccluded tube, p the fluid density, p the viscosity and p the pressure.
Mathematical modelling of an irregular arterial stenosis 3. GOVERNINGEQUATIONS
Blood flow through an axi-symmetric stenosis can be modelled in two dimensions using a cylindrical coordinate system. The dimensionless governing continuity and Navier-Stokes equations for this system are:
au, V, au,
a,+-+-“o, r
au, V, -jy
au, ap Y&=
--g
I +-
au, v, F+“’
az
au,
a3 I au --L+_-.L_‘+:
[ a9 ap I
Re
a5
ar r2
I
a%,
-=-jg+z a2
v
I
a221 au,
-g+;dr+j$[
,
(2)
a%, 1
’
1071
used (191 in the longitudinal direction and 19 radially), with a high concentration of elements near the wall and in the region of the stenosis (approximately 0.3 by 0.2 mm), decreasing with distance away from the stenosis. A system of algebraic equations is obtained from the mesh points and element connectivity by integrating the shape functions from the Galerkin weighted residual representation of the governing equations. The pressure degree of freedom has been eliminated using a penalty function formulation (parameter, E- 1 x 10T6) and is recovered after the velocities have been calculated. This has the effect of reducing the size of the system of equations. The system of equations is solved iteratively using up to three successive substitutions followed by quasi-Newton iterates with Boyden’s update. Generally, only five or six iterations were required for the solution to converge, leaving residuals of about 10e4 over the whole grid.
(3) where r and z are the dimensionless coordinates, scaled with respect to D,, with the z-axis located along the symmetry axis of the artery. No secondary or swirling flows have been allowed so the total velocity is defined by the dimensionless radial and axial components, u, and v,, scaled with respect to u,. These equations assume that the fluid is incompressible and Newtonian, which makes it reasonable to compare with the experimental results as water was the fluid used. Unfortunately, this theory does not lend itself to accurately modelling blood flow. The velocity boundary conditions on the arterial wall are the usual no slip conditions at r=R(z).
v,=vz=O
(4)
Flow upstream Ifrom the stenosis was assumed to have a velocity profile corresponding to Hagan-Poiseuille flow througb a long circular tube of constant crosssection v,=2
r2
( ) l--
,
R,2
v,=o,
(6)
(R, = D,,/Z). It was decided, after numerical study, to apply this boundary condBion eight arterial radii upstream from the beginning of the stenosis, thus avoiding interaction between the stenosis and the entering parabolic velocity profile. Downstream, the velocity is left free but conditions of &ro normal and tangential stress are applied. These ~conditions were applied eighty radii downstream from the beginning of the stenosis in order that fully developed flow be reestablished (Engtsman, 19% The govern+8 equations, being highly nonlinear, were solved nu+neri&ly using the FIDAP computational 5uid d ‘amice package, utilising the Gale&in bod (En&man, 1990,1982). A non5nite element F uniform grid of 13429linear quadrilateral elements was B?I 24:11-G
4. MODEL VALIDATION The model was validated by comparing modelled pressure drops with the experimental values given in the literature. Pressures are calculated by differencing the axial pressures at z=O.42 (0.13 cm) and z- 10.118 (3.12 cm) from the beginning of the stenosis. Note that these points are actually in the solid tubes attached at each end of the irregular section of the stenosis. A plot of the comparison is shown in Fig. 2. Good agreement is obtained for most of the range of Reynolds numbers considered. The difference is only of the order of a few per cent. Larger errors are observed with increasing Re and are due probably to the occurrence of 5ow instabilities or possible transition to turbulence (Deshpande et al., 1976), neither of which is accounted for in the model. These are regions where the 5ow patterns deviate from the smooth laminae assumed in the model. Over the physiological flow range (Re = 100400) the agreement is very good. This good agreement is the basis for using the model for a detailed study of the 5ow patterns obtained in the stenosis. Despite the fact that differences between the experiment and the model occur at higher Reynolds numbers, results are presented at these higher values for comparison with previous studies (Deshpande et al., 1976; For-rester and Young, 1970; Morgan and Young 1974).
5. SEPARATlONAND RFXIRCULATINCFLOW Points of flow separation are observed when the normal wall shear stress changes from positive to negative. (DetaiR of normal wall shear stress comput&ions arc given in Engleman, 1990.) The flow is said to reattach when the stress becomes positive again. A normalii wall shear stress is obtained by dividing calculated values by the value in the un-
P.R. iOHNSTON and D. KILPATIUCK
1072
II)
100 Reynolds Number
J 1000
Fig.2. Comparisonofmodelled and experimental pressure drops.
I
8,
I 11
7 -
I’
I' 11 ;I
6 -
Re-20 Re=lOOO ---- _
I’ I 1III !/I1
5 4 -
n
-1
’
0
I
2 4 Dimensionless
6 8 Axial Position
10
Fig.3.Wallshearstrewprolila along the stenosis at Reynolds numbers of 20 and 1000.
occluded tube section, far upstream from the stenosis (i.e. the value obtained from Poiseuillean flow). In the stenosis under consideration, there is no flow separation observed at low Reynolds numbers. At a Reynolds number of 600 a small region of recirculating flow develops between Z= 7.12 (2.19 cm) and z = 7.24 (2.23 cm) downstream from the beginning of the stenosis. Thii represents the separation Reynolds number, that is, the Reynolds number at which flow separation is first observed. When the Reynolds number is increased to 800, two recirculation zones are observed. The 6rst is a larger region in the same position as at Re=600, and the second is further downstream between z= 8.82 (272 cm) and z= 8.88 (2.74 cm). Increasing the Reynolds number further to 1000 produces another recirculation zone which is between the two at Re = 800, from z = 8.21(2.53 cm) to
z = 8.54 (2.63 cm) downstream from the beginning of the stenosis. It is interesting to locate the separation and reattachment points on the stenosis curve (refer to Fig. 1). The first region of flow reversal is in the large valley after the narrowest point of the stcnosia The second and third regions are in the next two valleys beyond the first. However, there is no recirculation in the fourth valley, near the exit of the stenosis. It is possible that a region may develop in this area at a higher Reynolds number. Figure 3 shows the wall shear stress vs axial position at Reynolds numbers of 20 and 1ooO. Note that the two curves are very similar in shape. Although the figure shows plots for only the extreme values of the Reynolds numbers considered, all curves for the intermediate Reynolds numbers lie between these two. The
Mathematical modelling of an irregular arterial stenosis
1073
where z,, is the half-length, 6 the width and z1 the centre of the stenosis. Here thii is used for comparison with results obtained from the irreguhu stenosis The centre of the smooth stenosis is taken at the narrowest point of the irregular one. Figure 4 shows a comparison of the smooth and irregular s&noses. A mesh containing 3629 linear quadrilateral elements, of similar size to those utiri above, was used to solve this model. 6. COMPARISONWllIi A SMOOTHSl’ENOSIg Comparing pressure drops over the two stenoses shows that the smooth stenosis predicts higher values Most model&g studies performed use a smooth than the irregular stenosis (Fig 5). The two curves are stenosis of the form parallel but the model of the irregular steno& gives much better agreement with the experimental results. A plot of the wall shear stress shows that the smooth (z)= stenosis predicts a much narrower range of values (Fig. 6). In the converging section of the stenosis the z>hJ, (7) irregular model predicts lower shear stresses, but at
points of separation and reattachment are clearly visible as the crossing points of the horizontal axis. Maximum shear stress occurs at the narrowest point in the stenosis. The expansion downstream of the narrowest point results in a rapid drop in shear stress in that region. In fact, the shear stress profile retlects very closely the outline of the stenosis.
l+u)s(y)), i:: loo0 (compared with 600 for the irregular model). It is also further away from the narrow point of stenosis than seen with the irregular stenosis. Another point of comparison is the axial velocity. Figure 7 shows that velocities through the irregular stenosis are less than those in the smooth one. However, it appears that the shapes of the velocity profiles are similar, implying that the tendency for flattened velocity proflles through the stenosis is the same in both models.
2
,4
7. MORE SEVEW SIXNOSES
To study the effects of a more severe stenosis, the original stenosis profile was modified to achieve a severity of 87% areal occlusion, while maintaining the same general shape. The results profile is shown in Fig 8, along with the original profile. A plot of wall shear stress is shown in Fig 9 for Reynolds numbers of 20 and 1000. These curves are smoother than those for the l& severe stenosis (Fig 3), indicating that the more severe stenosis tends to appear smoother in the diverging section. Also, the values are an order of magnitude larger than for the lesser constriction. There is only one region of flow recirculation which encompasses the three regions observed previously and it occurs for Reynolds numbers > 60. The presence of only one region is due to the
-8
6
Dimensionless Axial Position Fig. 6. Wall shear stress pro&s for irrqular and smooth models at Reynolds numbemof 20 and 1ooO.
‘\ ‘\ \
‘\ \
, I L
@
\\
\
\\
\
\
\
\
\
\
\
.35 .l .15 .2 .25 .3 Dimensionless Radial Position
0.05
.4
Mathematical modelling of an irregular arterial stenosis
1075
.s .4
.3
.2
Original Modified
.l
0
2
4
Dimensionless
Axial
6
Data Data
---_
8
10
Position
Fig. 8. Modified irregular stenosisand the original stenosis.
80
Re=lOOO ----
70 60 50 40 30 20 10 0 A”
0
2
4
Dimensionless
8
6 Axial
10
Position
Fig. 9. Wall shear stress profiles along the stenosis for the mcxlilkd irregular stenosis at Reynolds numbers of 20 and 1000.
Modified Smooth 0
’t 0
Original Data Approximation
----
1 I
2
4
Dimensionless
6 Axial
8 Position
Fig. 10. Smooth approximation to the mod&d original data
10
P.R.JOHNSON andD. KILPATRICK
1076
higher velocities achieved in the stenosis, allowing little opportunity for the flow to reattach to the walls even when there is a slight recontraction in the pro&. As previously stated, a smooth stenosis profile was developed to compare with this more severe stenosis. Figure 10 shows the two profiles. Rressure drop comparisons are shown in Fig 11. The two profiles are parallel again, as for the less severe stenosis, but the smooth profile model again overestimates the total pressure drop, this time to a greater extent. A plot of wall shear stress at Reynolds numbers of 20 and loo0 is shown in Fig. 12. Again the smooth profile model overestimates the stress in the converging section but underestimates it at the narrowest point. In the diverging section of the stenosis, the smooth model matches the irregular model very closely and, more importantly, gives a very good prediction of the point of flow separation. The separa-
tion Reynolds number for this model is < 60, which is an underestimate of the irregular model value. This is in contrast to the less severe situation, where the smooth model overestimated the separation Reynolds number. These observations are due probably to the similarity in stenosis pro&s in the diverging section. Finally, Fig. 13 shows axial velocity pro&s for the two severe stenoses. The smooth profile predicts higher velocities, as was observed in the less severe case. However, there is a difierence in the shape of the velocity profiles. The typical blunt velocity profile of flow through a stenosis is more exaggerated with the irregular profile. a
CONCLUSION
The mathematical model of flow through an irregular stenosis presented here shows excellent agreement
1000
Modified Irregular Model I Smooth Stenosis ---- ’
1.00:
.A”
~
10
100 Reynolds Number
Fig 11.Prcssuredropformodificdirrcgularsnd ~00th mod&.
no-20
g
g
ao-
lb-20
70.
ti
t 60?z 50 2 g 40. 3
(Irregulars
It,-1000 (Irregulars
It.4000
,swoth) @mothl
---1 I, ....'...I,'I ----
II
30.
fE 20%A 01 10 4
ca
-10 ' 0
2. 4 6 8 Dimensionless Axial Position
Fig. 12. Wall shear strcsa profrlcsfor modi!icd irregular nod smaoth mod& at Rcynold8 numbers of20 and1000.
Mathematical modelling of an irregular arterial stenosis
------------t
--__
1077
Irregular Model *_._\ Smooth Model ----
6
0 0.020.040.060.08 .l .12 .14 .16 .18 .2 Dimensionless Radial Position Fig. 13. Axial velocity profiles for the modified irregular and smooth models at z=6 and Re=200
with publish+ experimental data for the same stenosis. This made possible a detailed study of flow patterns inthkstenosis and the identification of three recirculation .xones in the diverging section of the stenosis. A comparison of results from the irregular model with the more common smooth model showed that the latter overlooked several fluid dynamical aspects of the irregular situation. These were overestimates for pressure droR, wall shear stress (except at the maximum) and separation Reynolds number. Also, only one recirculation zone was identified. These errors are due mainly to the simple shape assumed for the stenosis. When a more severe stenosis of similar shape was studied, only one recirculation zone was observed, due to the higher fluid velocity in the narrow section of the stenosis. A smooth stenosis of the same severity again overestimated the pressure drop and wall shear stress but underestimated the separation Reynolds number. It can be concluded from this study that it is necessary to have a good representation of the stenosis profile to make accurate predictions of the fluid flow in the stenosis. This is true especially when looking at recirculation xones. These are postulated to be prime areas for further deposition of atherosclerotic plaques which can lead to a worsening of the occlusion (Asakura and Karino, 1990). Methods based on fibre optic laser Doppler anemometry (Webber, 1982) and ultrasound (Kitney et al., 1990) are available to obtain estimates of the internal protIle of a s&nosed artery. These data could be used easily in a model, such as has been developed here, to give accurate predictions of the art&al fluid flow (assuming a gid artery with a Newtonian fluid) and therefore hi 4 light any potential areas of further deposition or extreme wall shear stress. Clearly, a threedimensional Buid flow model would give even better
insight, but this is not practical from a computational point of view at this stage. Other factors which may be important and should be studied are the inclusion of arterial wall motion and a non-Newtonian model for blood. Acknowledgement-PRJ was ,supported by the National Heart Foundation of Australia. REFRRRNCES
Asakura, T. and Karino, T. (1990) Flow patterns and spatial distribution of atherosclerotic lesions in human coronary arteries. Circulation Res. f&1054-1066. Back, L., Cho, Y, Crawford, D. and Culfel, R. (1984) Effect of mild atherosclerosis on 5ow resistance in a coronary artery casting of man. J. bionw?ck.Engng 1% 48-53. Deshpande, M., Giddens, D. and Mabon, R. (1976) Steady flow through modelkd vascular stenoses. J. Bimechics 9, 165-174. En&man, M. S. (1982) FIDAP-a fluid dynamics analysis &age. Ado. &I@ S&V. 4,163-166. _ _ Enakman M. S. 11990)FIDAP User’s Guide. Fluid Dvnamh International, Evanston, Illinois. Version 5.0 E&t. Forrester, J. H. and Young, D. F. (1970) Flow through a converging-diverging tube and its implications in occlusive vascular dVIZ theoretical development. J. Biomecka?tics3,297-305. Johnston, P. R. and Kilpattick. D. (1990) A mathematical model of paired arterial stenoses. In Computers in Curdiology (Edited by Ripky, K), pp. 229-232. IEEE Computer Society pnsS Chicago. Kit&y, R. I., Straughan, Kz Moura, L., Burr& C. and Rothman M. T. (1990) 3-d visuahsation for the studv of arterial disease and tissue chamcterimtion. In Com&ters in Cardiology(Edited by Ripley, K.), pp. 3-6. IEEE Computer Society Preq Chicago. Lee, J.-S. and Fung, Y.-C. (1970) Flow in looally constricted tubes at low Reynolds numbers. J. appl. MvcA 37,9-16. Morgan, B. and Young, D. (1974) An integral method for the analysis of flow in arterial stenoses. Bull. math. Biol. 36, 39-53. Webber, S. (1982) The haemodynamks of human coronary artery stenosis. M.Sc. thesis, University of Tasmania.
