Ind. Eng. Chem. Res. 1995,34, 3426-3436
3426
Convective-Diffusive Transport and Reaction in Arterial Stenoses Using Lubrication and Area-AveragingMethods? Marc Homer: Pedro Arce,* and Bruce R. Locke Department of Chemical Engineering, FAMUIFSU College of Engineering, Florida A&M University1 Florida State University, Tallahassee, Florida 32310
Atherosclerosis is a circulatory system disease characterized by the accumulation of cholesterol, calcium, complex carbohydrates, blood and cellular remnants, and fibrous tissue in the intimal layer of the artery. The deposition of these particles at the arterial wall may result in a stenosis, a decrease in the diameter of the artery. The focus of the present paper is to investigate the effect of the stenosis on the transport of the low-density lipoprotein, or LDL, from the bloodstream t o the arterial wall. This analysis models the hydrodynamics of the artery using lubrication theory. This method yields a Poiseuille type of velocity profile, with a modification due to the changing vessel diameter. Mass transfer of LDL particles to the wall of the artery is considered using the species continuity equation for dilute solutions. A pseudo-steady-state model is developed since the radius of the artery is assumed to decrease much more slowly with time than the rate required for the establishment of concentration profiles in the solution. Convective transport in the radial and axial directions couples the hydrodynamics t o the species balance equation. Diffusion in the axial and radial directions is also considered. The effect of varying convective, reactive, and geometrical parameters is analyzed. A first approximation, which assumed negligible radial convective contributions to mass transport, yielded quite different results compared to cases where this term was included. Thus, from the analysis, it appears that the convective transport in the radial direction is very important when considering mass transport in a channel of changing diameter, even when the relative change is small. 1. Introduction An arterial stenosis is the physical manifestation of atherosclerosis,a disease of large and medium arteries (Liepsch, 1990). This circulatory disease is the leading cause of death in the United States, as well as most other industrialized nations. Present atherosclerosis research is conducted in various disciplines including genetics, molecular biology, and engineering (Nerem, 1992). The possibility that the hydrodynamics present in the circulatory system may be a contributing factor in the growth and development of stenoses is just one area of engineering concern. The present research considers the connections between stenosis formation and the hydrodynamics present in the circulatory system. The focus of this investigation will be primarily on the effects of convective-diffusive transport of the low-density lipoprotein in blood. We believe that this study will be helpful in improving the understanding of the stenosis phenomenon and in determining “early” trends in the deposition of material on the arterial wall. The background for this investigation begins with cholesterol, a key component of all eukaryotic plasma membranes, which is essential to the growth and viability of cells in higher organisms (Stryer, 1988). Cholesterol biosynthesis in cells is one source of cellular cholesterol requirements. Biosynthesis is a complex process initiated by acetyl CoA in the cytoplasm of liver cells. Cholesterol can also be obtained from the diet in the form of the low-densitylipoprotein,or LDL, the chief transporter of cholesterol in blood (Stryer, 1988). The LDL particle is a macromolecule containing approximately 1500 cholesteryl esters. These molecules + Paper submitted by invitation on the occasion of the 35th Anniversary of Transport Phenomena by R. Bird, W. Stewart, and E. Lightfoot. * Present address: Department of Chemical Engineering,
Northwestern University, Evanston, IL 60208-3120.
0888-5885/95/2634-3426$09.0OlO
are surrounded by a hydrophilic envelope of phospholipid and unesterified cholesterol molecules. This coat also contains one copy of apoprotein B-100, a signal protein involved in the uptake of LDL into cells. The cell can use dietary cholesterol in membrane biosynthesis or store it for future use. Any deficiencies in the membrane biosynthesis can lead to decreased cellular uptake and possible increased concentrations of LDL in the circulatory system (Stryer, 1988). This can prove harmful, even fatal, to humans. Atherosclerosis is the primary example of a circulatory disease that may occur when the level of LDL in blood is too high. This disease manifests itself as an arterial stenosis, a progressive occlusion of the arterial conduit. The stenosis represents a focal point for accumulation of cholesterol and fat in the intimal layer of arterial wall tissue (Fairbairn et al., 1972),and can grow throughout life by continual deposition of material at the arterial wall. It has been suggested that the first stage in the development of atherosclerosis is the formation of fatty streaks (Hurst, 1990). These streaks consist of a collection of lipid-laden intimal cells which have two fates in the atherogenic process. They can either return to normal tissue, a reversible process, or they can develop into a fibrous plaque. The latter case is the first dangerous step in the atherosclerosis process. The fibrous plaque is composed of intimal lipid deposits which form a deep core of extracellular lipid and necrotic cell debris. This core, in turn, is covered by a fibromuscular cap (Glagov et al.,1990). The continued development of the fibrous plaque can lead to serious illness, possibly even death, in diseased persons through complete occlusion of the arterial conduit. Blood is the medium for transport in the circulatory system, and consideration must also be given t o its behavior in the arterial network. Whole blood is a suspension of erythrocytes, leukocytes, and platelets in plasma, and it is regarded as a colloidal dispersion
0 1995 American Chemical Society
Ind. Eng. Chem. Res., Vol. 34,No. 10, 1995 3427 because of the particulate nature of the red blood cells (Lightfoot, 1974). Plasma is a solution of salts and small proteins that behaves as a Newtonian fluid. Although blood is considered a colloidal dispersion from the physicochemicalpoint of view, its fluid behavior depends on the diameter of the flow conduit. This phenomenon was discovered by Fahraeus and Lindqvist (1931),who showed that the viscosity of blood is constant in tubes of diameter greater than about 0.3 mm. However, blood viscosity decreases as the diameter of the tube decreases below this value. This behavior of blood is known in literature as the ”Fahraeus-Lindqvist effect”. This response is due to a decrease in vessel lumen diameter which results in an increase in cell-wall interactions, causing a shift t o non-Newtonian behavior (Fahraeus and Lindqvist, 1931). The rheology of blood was f i s t explored by the French physician Poiseuille in the 1840s (Fahraeus and Lindqvist, 1931). Due to the experimental difficulties that arose when blood was used, Poiseuille further studied fluid behavior in glass capillary tubes using various other fluids, such as water. These early tests led to correlations between eMux time of fluid from a tube as a function of tube length, tube diameter, and pressure. k4 already mentioned, Fahraeus and Lindqvist showed that the behavior of blood in narrow capillary tubes changed from Newtonian to non-Newtonian as the tube diameter decreased below 0.3mm. This response was attributed to collisions between red blood cells and the channel walls, as well as t o the formation of rouleaux, aggregates of red blood cells, in the flow channel. These aggregates also contribute to the yield stress of blood. This yield stress was demonstrated by Merrill et al. (1965),in tubes 523 and 850 pm in diameter. The tube diameters used by Merrill et al., however, are far below the limits of the present investigation. After a pressure difference was applied across the end of each tube, a finite pressure difference was present across the ends of each tube when flow stopped. This pressure difference represents the yield stress of blood, and was shown to decrease with an increase in the diameter of the tube. Experiments where zero flow was approached from a flowing situation were performed because the settling of red blood cells from suspension eliminates the quantification of this yield stress from blood initially a t rest (Wayland, 1967). A 35 year study by DeBakey et al. (1985)demonstrated that atherosclerosis is a highly localized disease in the human circulatory system. The primary sites for the formation of stenoses are a t major branches of the aortic arch, visceral branches of the abdominal aorta, the coronary arterial bed, and the terminal abdominal aorta and its major branches. Fluid mechanical factors contributing to atherosclerosis a t arterial bifurcations, or branch points, include the development of regions of low and high shear, as well as the occurrence of separation points. It has been suggested that transport of solutes in areas of low shear is inhibited which may result in an increase in the accumulation of lipid in the arterial wall (Caro et al., 1970). The recirculation present in regions of flow separation may also contribute to an increase in residence time of thrombogenic material (Zarius et al., 1983). It has also been shown that regions of high shear, on the order of 400 dyn/cm2,in the canine aorta can cause endothelial injury, and thereby leave the arterial wall susceptible to particle entry (Fry,1968). These factors are also believed t o contribute to stenosis formation in arterial conduits. The fluid mechanical aspects present in an arterial stenosis have been studied using a number of ap-
proaches. Theoretical studies include that of Lee and Fung (19701,who investigated flow in a restricted channel using conformal mapping. This investigation was conducted for an axisymmetric stenosis with a 50% reduction in radius. Their results for a Newtonian fluid revealed an asymmetric velocity profile for finite Reynolds number. Misra and Kar (1989)explored blood flow through a stenosed blood vessel using a momentum integral approach. This study included the effect of slip a t the wall, a characteristic believed present under certain conditions of blood flow (Bennett, 1967). It was shown that this slip velocity increased the axial velocity through the stenosis. Young and Tsai (1973)investigated flow characteristics in axisymmetric and asymmetric stenosis models. They considered reductions in radius ranging from onethird to two-thirds that of the parent channel. Their investigation of separation points in the axisymmetric stenosis revealed that a critical Reynolds number of 190 was required for fluid separation in the stenosis of onethird reduction in radius. Stenoses with a two-thirds reduction in area had a critical Reynolds number of 10. Their analysis also showed that Poiseuille’s law gives an order of magnitude underestimation of the pressure drop. This difference is due to the increase in viscous loss that occurs when a constriction is present in the flow channel. The significance of this loss is important a t low Reynolds numbers where flow is very serious, and it is not accounted for when using the Bernoulli equation. Previous investigations concerning the solution of convective-diffusive equations in an arterial stenosis primarily focused on the supply of nutrients, such as oxygen, to the arterial wall. Schneiderman et al. (1979) analyzed the transport of oxygen to determine whether or not tissues were affected by the presence of a stenosis in the arterial conduit. Regions of flow separation were found t o have decreased levels of mass transport, while vortex regions actually showed some enhancement. The hydrodynamic solution of Morgan and Young (1974)was used t o specify the radial and axial velocity components. This model used an integral approach to yield an approximate analytical solution of the Navier-Stokes equations. Nazermi and Kleinstreuser (1989)analyzed particle trajectories in an arterial bifurcation by looking for hemodynamic patterns that may be responsible for increased wall deposition. The results of their numerical solution showed that particles entering near a wall will most likely impact on the wall in regions typically viewed as susceptible to atherosclerosis. It also showed lateral migration of particles due t o diffision to have a negligible effect. It appears that the application of mass transfer results using locally derived fluid mechanical models through an arterial stenosis have not been reported in the literature, and will be the focus of this paper. A geometrical model is used in the present work to describe the stenosis in the arterial conduit. This model is based on two inverted cylindrical cones whose diameter is a parameter in the analysis. It is believed that the hydrodynamics present in the artery may play a key role in the growth of arterial stenoses, and it is the effect of the hydrodynamics, in conjunction with the convective-diffusion transport, that will be investigated here. This will begin with the development of a differential model capable of describing the velocity profiles in the stenosis domain. This model will be derived using lubrication theory as the starting point. Then, the mass transport process will be described using the molar form
3428 Ind. Eng. Chem. Res., Vol. 34,No. 10,1995
of stenosis. H(z) is a discontinuous function a t the origin, z = 0. This requires treating the stenosis as two separate sections, one to the left (i.e., the convergent section) and one to the right (i.e., the divergent section) of the origin. These two regions are coupled a t the origin, their point of intersection. The mass continuity equation for an incompressible, Newtonian fluid under steady-state and creeping flow in cylindrical coordinates is (Bird et al., 1960)
+ L l
+.L+
(2) The Navier-Stokes equation of motion in the rdirection for the system described above is (Bird et al., 1960) IbJ
Figure 1. Axisymmetric stenosis system.
of the species continuity equation for dilute solutions. The convective-diffusive mass transport equation is sequentially coupled to the equation of motion. An effective transport model describing the concentration profile of LDL particles in the stenosis domain will be derived using an area-averaging approach as developed by Slattery (1981)and Whitaker (1985). The resulting averaged model is simplified further in order to obtain a first approximation to the concentration pmiile through the stenosis. Emphasis will be on the early stages of stenosis development, and the hydrodynamic model will assume the behavior of the blood in the limit of Newtonian fluids. The effect of stenosis geometry, the wall flux, and the radial velocity component on the concentration profile will be analyzed through the solution of the species continuity equation. 2. Formulation of the Hemodynamic Model
The hemodynamic model for the stenosis system, as predicted by lubrication theory, is derived using the procedure given by Denn (1980). The geometry and results of this analysis are, however, new applications of lubrication theory. The system studied here is an axisymmetric stenosis, as shown in Figure 1. This figure depicts a schematic of the model used as a n idealization of true stenoses found in the human circulatory system. Note that a scaling factor, x , is used to denote the radius of the flow channel a t the stenosis throat as a fraction of the radius at the stenosis entrance. The system is modeled using cylindrical coordinates, where the assumptions for the fluid mechanical model are that the 8 component of the velocity is zero, well-developed, laminar flow is present, there is no slip at the walls of the artery,there is incompressible flow, and there is fluid with constant viscosity. Satisfaction of the qong-flow channel” and “nearly parallel” assumptions is also required if a lubrication analysis is to be utilized. These state that R/L and (R - xRKH(z) must each be much less than unity, respectively. All the symbols have been defined in the Nomenclature section. The width of the artery, Hk), is determined using the following equation:
Equation 1 specifies each side of the stenoses (i.e., convergent and divergent sections) in the proposed linear form. Although this function is relatively simple, it represents fairly well the early stages in the formation
and the z-component of the Navier-Stokes equation is
Equations 2-4 are made dimensionless by scaling each variable to its characteristic value. These scaling values are R and L in the radial and axial directions, respectively, x is the characteristic pressure, and V and U are the characteristic velocities in the radial and axial directions,respectively. While the characteristic lengths in each direction are known, the characteristic velocity and pressure must be determined by using these values. The scaled continuity equation in conjunction with the z-component of the equation of motion (Denn, 1980) yields V = (R/L)U and x = pUL/R2,respectively. This is the appropriate form for the characteristic pressure associated with viscous flows. Substituting for the characteristic pressure and for the axial characteristic velocity into the r-component of the dimensionless equation of motion reveals that the pressure gradient in the r-direction is negligible with respect to the other terms in the equation. This conclusion is baaed on the ’long-flow channel” condition. This also implies that the approximation used in deriving the velocity profiles is more accurate for stenoses in the early stages of development, i.e., the kearly parallel” situation (Homer, 1994). The z-component of the dimensional equation of motion has the following form:
(5) where the pressure is a function ofz only. The boundmy conditions for eq 5 are B.C.l:
%HC*,
B.C.2
q,=o
=O
-
a bounded function
(6)
(7)
Equation 6 represents the no-slip condition at the wall, and eq 7 is a statement of the continuity of momentum flux a t the center of the arterial conduit. The velocity in the z-direction is found by integrating eq 5 and applying the above two boundary conditions, to obtain
Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 3429
-0.005 "r
,,,z" -0.01
---0 0
Mid-Point Between Entranceand Throat Stenosis Throat I
I
I
0.2
0.4
0.6
.\
.--. ---
Near the Stenosis Eninnce
'* 0.8
Mid-Point Between Entmnce w d Throat Stenosis Throat
1
-0.015
0
1
I
1
1
0.2
0.4
0.6
0.8
1
A
Figure 2. Typical axial velocity profiles in the stenosis domain (R= 0.3 cm, r R = 0.27 cm, L = 1.0 cm).
Equation 8 is similar to the Hagen-Poiseuille result (Bird et al., 1960). The difference lies in the fact that the axial velocity is now a function of the changing channel radius, i.e., H(z). The ratio of the axial velocity t o the maximum velocity is
r
Figure 3. Typical radial velocity profiles in the stenosis domain (R= 0.3 cm, XR = 0.27 cm, L = 1.0 cm).
in the divergent section, where the walls of the artery return to the original width (Horner, 1994). The pressure gradient in the stenosis can be determined implicitly as a function of the volumetric flow rate. This is done by first determining the volumetric flow rate through integration of the velocity profile by (11)
Typical numerical values for the volumetric flow rate can be found in the literature, once the region of interest in the body is identified. The pressure gradient can now be determined by Equation 9 relates the axial velocity to the axial position in the convergent-divergent sections of the stenosis. This analytical expression will prove valuable in the analysis of the convective-diffusive transport in a stenosis domain. This study will be conducted in a later section of this paper. Typical profiles given by eq 9 are shown graphically in Figure 2. This figure shows that the dimensionless axial velocity changes with position down a stenosis where R is 0.30 cm, x is 0.9, XR is 0.27 cm, and L is 1.0 cm. The axial velocity increases down the convergent stenosis section for the case analyzed. This is expected since the area available for flow decreases as percent stenosis increases. A relationship between the axial velocity and percent stenosis also can be deduced from the figure in case it is needed (Horner, 1994). Percent stenosis is defined as the percent reduction in the crosssectional area of the open artery at the stenosis throat to that of the unaffected artery, and is equal t o (1 - x ) x 100 here. Analysis performed on the velocity profile for different percent stenosis reveals that the axial velocity increases as percent stenosis increases (Horner, 1994). This also is attributed to the changing crosssectional area of the flow channel. The average axial velocity will be required for the area-averaging analysis in the mass transport section. The average velocity is found by integrating the axial velocity over the cross-sectionalarea of the stenosis, and then dividing by that area (see eq 8). This procedure yields the following average velocity
Note that the averaging operation removes the rdependence from the axial velocity. Equation 10 indicates that the axial velocity increases as the convergent section of the stenosis is traversed, and then decreases
(12)
Finally, the velocity in the r-direction can be found using the mass continuity equation, eq 2, in conjunction with the axial velocity, eq 8. Substituting for the pressure gradient using eq 12 allows elimination of the z-dependence on pressure in this analysis. Integration of the resulting continuity equation yields,
Equation 13 predicts that the velocity in the r-direction is zero if the channel has straight walls. It also predicts that the radial velocity a t the arterial walls is zero, in accordance with the no-slip condition used. Equation 13 may be put into more familiar form by substituting for the volumetric flow rate using eq 12 to obtain
Typical plots illustrating the behavior of eq 14 are presented in Figure 3 for x equal t o 0.9. The velocity profiles are shown in nondimensional form by using the ~ and , the ~radial ~ coordinate P. maximum velocity u This figure shows a parabolic type profile where the radial velocity is zero at the wall, in accordance with the no-slip condition used in this model. The vertex of the parabola moves toward (radial) positions closer t o the wall when the axial position moves from a point near the stenosis entrance to the midpoint and at the stenosis throat, respectively. The position of the vertex covers a relatively small range, however, located in the range 0.5 s rlH I0.6. There is also no "inversion" between the branches of the cases mentioned above for the values of rlH before the vertex. The point of inversion is
3430 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 located around r/H = 0.7. The maximum of the velocity component given by eq 14 shows the position of the maximum velocity t o occur at a radius of [H%)/311’2. Comparison between profiles for different percent stenoses shows an increase (Horner, 1994) in the radial velocity as percent stenosis increases. Physiological Considerations. Certain limitations are present in the hydrodynamic model derived above due to the use of the lubrication approximation to simplifj. the Navier-Stokes equations. Of primary concern is fulfillment of the laminar flow requirement, and a secondary aspect is that the reduction in channel radius must not exceed 20% (Denn, 1980). It is generally accepted that atherosclerosis is a localized disease; i.e., certain areas of the human circulatory system are more a t risk than others. The study by DeBakey et al. (1985) revealed atherosclerotic lesions to be most likely found in one of the following four regions: the branches of the coronary arterial bed, the major branches of the aortic arch, the visceral branches of the abdominal aorta, and the terminal abdominal aorta and its major branches. The first three categories represent buildup in the area of an aortic bifurcation. Here it is indeed believed that hydrodynamics has a direct effect on intimal thickening of the arterial wall (Zarins et al., 1983). Analysis of blood flow in a bifurcation has revealed the presence of regions of low as well as high shear stress. Regions of low shear typical occur downstream from the stenosis throat. Such factors as flow separation and recirculation appear to contribute to lipid deposition in these areas. At high shear rates, it has been suggested that mass transfer of lipid to the wall increases (Nazemi and Kleinstreuer, 1989). It is the latter category of the DeBakey study where the proposed model of this paper is more appropriate. More specifically, the femoral artery is, perhaps, the section of the circulatory system where the model is most applicable. The flow conditions expected in the femoral artery must now be identified. One way to obtain an idea about the flow regime is to calculate the Reynolds number, and use this value to justify the laminar flow assumption made in the hydrodynamic analysis. The Reynolds number for the system is defined as N R=~cud/ p. This definition is based on the tube diameter, a standard procedure when considering flow in a pipe or tube. The inner diameter of the femoral artery is approximately 6 mm (Dawson, 19911, and the mean velocity there is on the order of 5 cm/s (Dawson, 1991). The apparent viscosity of blood is 0.040 g/(cm s) in the aorta, and the density of blood is 1.05 g/cm3 at normal hematocrit (Pedley, 1980). Using these values, the Reynolds number is found t o be approximately 75 in the femoral artery, which clearly indicates that the flow regime is laminar. A Reynolds number of 75 is also far below the critical Reynolds number for flow separation of 190, which was stated by Young and Tsai (1973) for an axisymmetric stenosis with one-third reduction in area. Under such conditions,the flow can be considered primarily viscous, and with inertial terms of negligible value (Perry and Green, 1984). The next consideration is the physical characterization of the LDL particle. As stated in the introduction section, the LDL particle is made up of about 1500 cholesteryl esters coated with a layer of phospholipid, unesterified cholesterol, and one copy of apoprotein B-100. The LDL has a diameter of 22 nm and a molecular weight of 3 million g/mol. The concentration moVcm3 in the average of LDL in blood is 2.27 x human (Stryer, 1988). The small diameter and low
concentration of the LDL particle suggest a continuum approach to the mass transfer of LDL particles in blood can be utilized. 3. Mass Transport of LDL in the Stenosis Domain The molar form of the species continuity equation will be used to describe the mass transfer of LDL within the stenosis domain. The solution to this equation is derived using the area-averaging approach developed by Slattery (1980), Whitaker (1985) and others. The advantage of this approach is the elimination of the radial dependence from the species continuity equation, through identification of effective transport coefficients. These effective transport parameters help in understanding the convective-diffisive transport in the stenosis region. The consumption of LDL particles at the arterial wall will be modeled using zero- and firstorder reaction rates, which lead to a Robin boundary condition at the surface of the arterial conduit. The zero-order condition assumes a constant value for the flux of LDL particles a t the arterial wall, whereas the first-order boundary condition is a function of the local LDL particle concentration. The concentration field of species in a homogeneous fluid, neglecting bulk reactions, is found using the molar form of species continuity equation (Bird et al., 1960)
at
+ V.(vC)- DV2C = 0
(15)
where convection and diffision were considered the two sources for flux of the LDL particles in solution. The system is modeled by using cylindrical coordinates under the assumptions that the blood can be assumed to be a homogeneous solution, LDL concentration is dilute, there is no angular dependence on concentration, and the angular component of the velocity is zero. The model also assumes constant density and diffisivity, and the following analysis is valid only for steady state. In order to apply the species continuity equation, it is assumed that the particles move at the local fluid velocity. The species continuity equation (15) for the case of cylindrical coordinates reduces to
under the assumptions stated above. Also, the dilute solution assumption allows the mass average velocities developed in the hydrodynamic section to be utilized in place of the molar average velocities required in eq 16. This equation is recast in nondimensional form using the same dimensionless variables as in the hydrodynamic analysis performed in the previous section. The only additional definition is the dimensionless concentration which is defined as the ratio of the concentration a t any point to that at the inlet. The resulting dimensionless equation can be simplified further by substituting for the characteristic velocity, V, and by defining *
UR and & = -R L
Pe =-D
17)
to give
18)
Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 3431 The convective terms in eq 18 carry much more weight than the diffusive terms, since the coefficient of the nondimensional diffusive terms is on the order of under the physiological conditions considered here. Thus the scaling analysis predicts the importance of convection in the transport or LDL within the stenosis domain. Three cases will be considered using this approach. Cases 1 and 2 assume negligible effects of the radial velocity relative to the axial velocity component on the mass transport of LDL within the stenosis domain. Zeroand first-order reactions at the wall account for deposition. For the zero-order reaction, the development of an averaged equation will give effective transport parameters, namely diffusion and velocity terms, that can be used to get qualitative and quantitative information on the transport characteristics of the system without a complete solution t o the original partial differential equation. F'urthermore, the averaging analysis is used to simplify the convective-diffisive equation describing the system. Finally, case 3 considers the effect of radial transport through the inclusion of the radial convective term in eq 16. This step-by-step procedure is used here so that the impact of each term on the concentration profile can be examined. 3.1. Area-AveragingSolution to the Differential Model. The first step in the area-averaging approach is to define the average by using (Bird et al., 1960)
where f is an arbitrary function. Therefore, the rdependence on the concentration will be eliminated in the averaged quantity 0. The first two cases to be examined assume that the axial convective transport is dominant with respect to the radial convective contribution and that the flux of the LDL a t the wall is either constant or a function of the local fluid concentration. The effect of the radial convective transport will be studied in a later section. The convective-diffusive equation to be averaged in this case is given by
a s) ac + g a% ]
ac
u - = D[--dr 1 ,a2
condition typically used in chemical reactor design, which sets the flux equal to zero at the stenosis exit. Equation 24 can be regarded as the point where the artery is no longer susceptible to increased rates of LDL uptake. Taking the average of eq 20
The first term in eq 25 is expanded using Gray's decomposition (Gray, 1975)applied to the velocity and concentration fields u, = (u,)
+ 6,
(26)
likewise,
C=(C)+C
(27)
where the tilde denotes the difference between the average and the point value of the function. These terms are commonly referred to as the "deviation" between the point value and the averaged quantity. Using Gray's decomposition and the centerline boundary condition yields
Note that the average of the axial diffusion term is the average, since there is no radial dependence of the average. To determine the concentration profile for the average concentration in the stenosis domain, the deviation in the velocity must next be identified. Equations 8 and 10 in conjunction with eq 26 yield the deviation field for the velocity
[
€z)l
8, = (u,) 1 - -
(29)
Next, the deviation field for the concentration can be determined by first subtracting the averaged balance from the point balance to give
(20)
a2C
The following boundary conditions apply: B.C.l:
(21)
=O = kCn
B.C.2: B.C.3: B.C.4:
(22) (23)
- D zacl
Neglecting the gradient of the deviation in the axial direction gives
=O z=L
The two conditions required for the r-direction are represented by eqs 21 and 22. These are the continuity of LDL flux a t the center of the channel and at the wall, respectively. Notice that this model can account for lipid deposition through a specified constant flux (i.e., n = 01,or as a function of the concentration of LDL in the blood stream, n = 1. The two axial conditions are eq 23, the inlet LDL concentration, and eq 24, a
A constitutive equation that relates the deviation field of the concentration to the average concentration is now required to obtain closure (Paine et al., 1983). The following constitutive equation analogous to that proposed by Paine et al. (1983)can be used to describe the deviation field for the concentration: (32) where n is the order of the reaction a t the wall. This equation applies only to the two cases of zero-order and first-order reaction at the wall. The functionsflrp) and
3432 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
g(r,z)must be determined using the boundary conditions and the deviation equation. 3.2. Analytical Solution for Zero-Order Reaction at the Wall (n = 0). For the case where the consumption of LDL a t the wall is independent of the local fluid concentration, the function flr,z) can be determined by substitution of eq 32 into eq 31 and integrating twice. Noting that the derivative of flr,z) with respect to r must be finite at the centerline, one finds
an effective diffisivity to relate the contribution of radial diffision to mass transfer
where De, can be identified by comparison of eqs 39 and 37 as (40)
where it was also noted that the integration may be carried out directly, since dfldr is a given function a t the wall (i.e., r = H(z)). The average value of flr,z) is zero. This is found by noting that the average of the deviation in concentration is zero, and through substitution of this result in eq 32. The integration constant, ,8, can now be found by averaging eq 33 and setting it equal to zero. This results in the following equation for
This equation is the converging-diverging channel wall equivalent t o the Taylor-Aris result (Aris, 1956). The modification is that the Peclet number, and thus the effective diffisivity, is now a function of axial position in the stenosis. Calculation of DedD shows it t o be on the order of 1013under typical physiological conditions in the femoral artery. This allows elimination of the first term of eq 40, resulting in the following approximate equation for the average concentrationthrough the stenosis:
Ar,z):
Equation 41 is written in a nondimensional form by scaling z t o the length of one stenosis section, L , and C to inlet concentration, CO,i.e., 2 = z/Land = C/Co
e
The derivative of A r j ) evaluated at the wall, is zero. Therefore, the final equation for the deviation field of the concentration is given by where
where H ( z ) is the given function of z , and the average velocity is given by eq 10. An equation similar to eq 35 was first obtained by Taylor (1953) using a completely different approach and for the case of straight circular tubes. The last term on the right hand side of eq 28 can now be evaluated using eqs 29 and 35: (36)
Substitution of this result into eq 31 and use of eq 10 yields the following equation for the average molar species concentration:
The Damkohler number, Da, is the ratio of reactive transport in the system to the transport by diffusion. This group has different forms depending on the reaction order, n. For example the group will not include the inlet concentration if a first-order reaction were assumed at the wall ( n = 1). The geometrical factor Ge given by eq 43 is the ratio between the position of the stenoses wall, H(z), and the length of one cylindrical section of the stenoses domain, L. An analytical solution of eq 42 can be found by rewriting such an equation as
(44) By integrating twice the following equation for the concentration profile of the low-density lipoprotein in the stenosis domain results:
where Pe is a modified Peclet number defined by (38)
and a total derivative is now used since the average concentration is a function of z only. Equation 37 can be recast in the following form using
where the plus sign refers to the convergent section, the minus sign applies to the divergent section, and a and ,8 are integration constants.
Ind. Eng. Chem. Res., Vol. 34, No. 10,1995 3433 100
-
,
1
Da&.7E+3
1
60 : :1 70 I
x
x
I
I
1
I
0
0.5
1
Ped.OE+7 Pe=S.ZE+I Pe=6.7 E+l
50
0.2
0
0.6
0.4
1
-I
z^ Figure 4. Axial conversion profiles showing the effect of the Damkohler number (zero-order-reactioncase, Pe = 5.2 x lo7).
Equation 45 is subject t o the following boundary conditions in each region:
B.C.l:
(46)
B.C.2:
(47)
B.C.3:
d(ol
i Figure 6. Axial conversion profiles showing the effect of the Peclet number (zero-order-reactioncase, Da = 3.3 x lo3).
---
40 Y
I
1
0.5
1
ov -1
where I refers to the convergent section and I1 refers to the divergent section. Equation 46 specifies the inlet LDL concentration, eqs 47 and 48 represent the continuity of concentration and flux where regions I and I1 meet, and eq 48 sets the flux at the exit equal to zero. This last condition is an approximation in that the flux is some constant, k, in the stenosis domain, but then set equal to zero by this condition a t the stenosis exit. The above system of equations is solved simultaneously to determine the four constants of integration below:
1 €2 - (a a' = -Du L 12
+
Ge -
PI1=$
I
xR.0.30crn xR=Ll.21 crn
d(@"
(49)
a' + 1 - -Pe 48
-0.5
--Du-[ 1 Pe Ge R 12 48
(50)
---
24
]
- (d)
(53)
The concentration profile through the stenosis is now completely specified. Graphical illustrations for the average concentration profile are presented in Figures 4-6. The ordinate of each graph represents the uptake of LDL particles at the wall in terms of percent conversion of LDL particles in the blood stream. The conversion, x, is in terms of the percentage of LDL particles entering the stenosis. The abscissa is the dimensionless axial position. Each plot shows a steady, almost linear, increase in LDL conversion down the convergent section of the stenosis.
I -0.5
0
e Figure 6. Wall axial conversion profiles showing the effect of the stenosis extent (xR)(zero-order-reaction case, Pe = 5.2 x lo7,Da = 3.3 x 103).
However, the slope of the concentration profile decreases as the divergent section is encountered, with an almost flat profile a t the exit. The shape of the profiles can be interpreted as a combination of using the constant wall flux within the stenosis domain and the flux being zero a t the exit. This result seems t o indicate that the hydrodynamics present in the convergent section increases the rate of LDL uptake, while that present in the divergent section tends to reduce deposition. Figure 4 depicts a direct relationship between the magnitude of the Damkohler number and conversion, x. Thus removal of LDL from the bloodstream increases as the reaction rate increases relative to diffusion. Conversely, Figure 5 illustrates the inverse relationship between the Peclet number and conversion. Here, the conversion is limited by the time available for LDL particles to diffuse to the wall, leading to decreases in conversion as convection increases. Finally, Figure 6 shows the direct relationship between the extent of stenosis and LDL conversion. In this case, conversion decreases when the area available for flow decreases. This is consistent with the case of increasing Peclet since convection through the stenosis system increases when the area available for flow decreases. 3.3. Solution to First-Order Reaction Problem (n= 1). The gradient in the average species concentration a t the wall is now proposed to be a function of the species concentration a t each axial position. The constitutive equation describing the deviation field of the LDL particle also includes the local average species concentration, i.e., n = 1 in eq 32. The procedure outlined in section 3.2 yields the same function for fir&),
3434 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
while g(r,z) is now found to be
of the average of this field (both in the r-direction) are assumed t o be of the same order, i.e.,
Substitution for flr,z)andg(r,z) in eq 32 yields a new deviation field for the first-order case. Substitution of the new deviation field yields an ordinary differential equation describing the average species concentration for first-order reaction at the wall
This approximation is based on the assumption that the
where it was assumed that the Damkohler number, a reaction term described below, is small. Equation 55 case in nondimensional forms gives
The dimensionless groups Pe and Ge utilized in eq 56 were defined in section 3.2. The Damkohler number for first order reactions is given by Da = KLID. Equation 56 was solved numerically using the DBVPFD subroutine from the IMSL library, which solves a system of first-order differential equations by using the finite-difference method. Due to the flux of LDL particles at the wall of the channel, two boundary conditions are required for each equation. Here, eqs 56, 46, and 49 are reduced to a system of two differential equations and two boundary conditions. The continuity of concentration and continuity of flux boundary conditions are eliminated here through the use of the absolute value function in the computational scheme. The average concentration profile is also calculated aRer solving the equation for the concentration profile. Numerical results for the average concentration profiles were investigated graphically for the same values of Da, Pe, and percent stenoses used in section 3.2. These figures have the same qualitative slope as those of section 3.2, with numerical differences as described below. The reaction order affects the conversion values in the region zlL 0 (i.e., the entrance region) where the profiles show steeper slopes. Also, the conversion levels are less than for the zero-order-reaction case for different Da numbers. For example, the wall conversion, for the same values of Pe and Da of Figure 9 a t (zlL) = 1, is less than 30% as compared to over 40% for the case of zero order (Horner, 1994). 3.4. Mass Transport Including the Radial Velocity and First-Order Reaction. This analysis begins with the addition of the radial convective term to eq 28. Averaging this term and using Gray's decomposition yields
or (u,)
aC
- + 8, ar
aC
= (u,) ar
The addition of eq 58 to eq 31 yields the final deviation equation for this case. The procedure outline in section 3.2 is again followed. h a first approximation, the magnitude of the gradient of the deviation field and that
small radii found in the femoral artery will result in very small concentration gradients in the radial direction. Therefore, the average of the radial gradient of the deviation field is a good approximation of the gradient of the deviation. Use of eq 59 allows the elimination of the radial convective terms in the deviation field. This yields functions flr,z) and g(r,z) to be the same as in the previous case. Therefore, the constitutive equation describing the deviation field for the concentration is the same as that of the previous section. The average convective-diffusive equation, however, is modified by the two terms on the right hand side of eq 57. Evaluation of the additional terms, and the previously mentioned assumptions of small Damkohler and large Peclet, yields the following ordinary differential equation for the average concentration profile in the stenosis
It is seen here that the derivative of H ( z ) replaces the Peclet number group of eq 56. This implies that the effect of the stenosis slope on the flow field is much greater than the convective effects present in the system when the radial velocity is included. Equation 60 is recast in dimensionless form using the scaling factors proposed in section 3.3:
(61) The dimensionless variables and groups utilized in eq 61 were defined previously in section 3.1 and in section 3.2. The first step in the solution of eq 61 is transformation into two first-order differential equations as in the zero-order case. The convective-diffusive model is fully specified through utilization of the same boundary conditions as the previous problem. The solution is found using the DBVPFD subroutine in conjunction with the boundary conditions eqs 46 and 49 and the transformed differential equation for this section, eq 61. Numerical results for the average concentration profiles are presented graphically for varying Damkohler number, Peclet number, and percent stenosis in Figures 7,8, and 9, respectively. These figures show a considerable departure from the profiles shown in Figures 4-6 for the convergent region of the stenosis domain. The qualitative shape for the divergent region of the stenosis domain, however, follows the same trend shown in Figures 4-6. In particular, Figures 7-9 show reduced LDL uptake in the divergent section where the radial velocity is toward the center of the channel, and large increases in the divergent section where the radial velocity is toward the wall. Also seen is the fact that the trend in changing the percent stenosis is reversed in Figure 9, when compared, for example, to Figure 6. The LDL conversion increases as the stenosis slope, and thus the magnitude of the radial velocity, increases. The
Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995 3435 20
I
-
x
I
I
1
Da=6.7€*6
I
-1
-0.5
1
I
I
0
0.5
1
e Figure 7. Axial conversion profile showing the effect of the Damkohler number where the radial convective transport is included (first-order reaction case, Pe = 5.2 x lo7). 20
-
15
-I.-
I
I
I
PeJ.OE+7 Pe=S.ZE+7 Ped,7E+l
-
10 5 0
-I
-0.5
0
0.5
1
2 Figure 8. Axial conversion profiles for the case of Figure 7 showing the effect of the Pe number (Da= 3.3 x lo6). 20
I--
x
15
I
1
I
xR.0.24cm xR=O.ncm xRAMO cm
10
5 0 -5
L -1
I
-0.5
1
8
1
0
0.5
1
2 Figure 9. Wall axial conversion profiles for the case of the Figure 7, showing the effect of the stenosis extent (xR) (Pe = 5.2 x lo', Da = 3.3 x lo6).
differences observed here seem to support the importance of the radial velocity component in the transport of LDL encountering arterial stenoses. 4. Summary and Concluding Remarks
The results of this paper are useful in gaining a deeper understanding of the transport behavior present in the early stages of an arterial stenosis. The analysis features two main aspects in the convective-diffusive transport of such a problem. The first is centered on the description of the hydrodynamics of blood flow within the stenosis domain. The second aspect is
related to the prediction of the concentration fields through the stenosis. The first objective is accomplished using the lubrication approximation in a cylindrical (convergent-divergent) system. This method yields a relatively simple solution for the radial and axial velocities, which can be used to describe the transport of LDL particles in the stenosis system. The second objective is accomplished through formulation of a convective-diffusive transport model based on the species continuity equation coupled with chemical reaction at the wall of the arterial conduit. This microscopic model is then averaged to obtain an effective transport equation in the axial direction. The solution t o the resulting model depends upon the differential equation, where analytical and numerical solutions were developed and concentration profiles illustrated for a variety of parameter values. To the best of our knowledge, this research represents the first effort in applying the lubrication approximation coupled with averaging methods to the analysis of arterial stenosis problems. Lubrication theory predicts a modified Poiseuille flow axial velocity profile through a mild arterial stenosis. The modification is carried through as an axial velocity that is now a function of the changing wall radius H(z). The average axial velocity is inversely proportional to the square of the channel radius. This velocity is symmetrical with respect to each region, increasing down into the convergent section, and decreasing out of the divergent section. The model predicts values for the radial velocity that are approximately 1/100of the maximum axial centerline velocity. This is expected when the depth of the stenosis into the artery is small. The assumptions of slow flow in the system supports the symmetry. Turbulent effects such as flow separation and circulation may create irreversible, asymmetric conditions in the system but are assumed absent here mainly based on the assumption that the model describes a relatively mild stenosis situation. Furthermore, physiological parameters for this system dictate that the average Reynolds number for this system is approximately 75, which implies that a laminar flow regime is assured. The analysis shows how the presence of stenoses in the femoral artery affects the relative importance of the radial and axial velocities of the system. For example, the decreased cross-sectional area acts to increase the axial convection of blood through the stenosis, thus reducing residence time of LDL particles in the system. A parametric study shows how the convective transport affects the conversion of LDL by increases in the Peclet number which decreases the time available for reaction for the various values of the Damkohler number. Furthermore, the analysis of section 3.4 reveals the importance of the radial velocity component in the convective transport. The final result is in agreement with the scaling analysis presented in the beginning of the mass transport section. Within the assumptions, the convective-diffusive transport model predicts trends in the behavior of the concentration field and also indicates the possible role of the transport in the deposition of the particles on the arterial wall. Further investigation is needed to compare, for example, the effect of the convective transport beyond the limit of the lubrication approximation. The process of deposition, diffusion, and reaction at the arterial wall needs improved characterization. In this analysis, such a process has been modeled using a zeroand first-order (global)rate with a wide range of values for the Damkohler number. Some caution must be used in the analysis of the quantitative results of this paper
3436 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
since experimental data about the kinetics of the deposition process are not currently available. Also, further analysis is needed t o check the validity of some of the approximations in the process of obtaining effective transport coefficients in the stenosis domain. The effect of diffusion in the deposited layer (at the arterial wall) should be coupled with the transport in the “fluid” and the overall system studied. Finally, the motion of red blood cells and their effect on the mechanism of transport in the stenosis domain should be investigated. A preliminary analysis of this aspect is currently in progress.
Acknowledgment The work described in this contribution is based on a Senior Honors Thesis in the Major (Chemical Engineering) by Marc Horner. Portions of the work have been presented a t the Annual Meeting of the American Chemical Society, Florida Section, Orlando, FL, May 1994,and the Annual Meeting of the Division of Fluid Dynamics of the American Physics Society,Atlanta, GA, November, 1994. One of us (P.A.) was partially supported by the Council for Research and Creativity, Florida State University, during the research conducted for this paper. Final revisions of the draft were typed by Heidy Booeshaghi.
Nomenclature C = molar concentration of LDL, mol/cm3 D = diffusivity of LDL, cm2/s L = length of one section of the stenosis, cm R = channel width at the stenosis inlet, cm U = characteristic velocity in the axial direction, c d s V = characteristic velocity in the radial direction, c d s x = fractional radius of open artery at the stenosis throat, dimensionless xR = channel width at the stenosis throat, cm p = fluid viscosity, g/(cm s) JC = characteristic pressure, dimensionless x = reactant conversion, % = deviation field in Gray’s decomposition 1 1 = absolute value ( ) = averaged quantity = dimensionless quantity
-
A
Literature Cited A r i s , R. On the Dispersion of a Solute in a Fluid Flowing Through a Tube. Proc. R . Soc. London, A 1966,235, 67.
Bennett, L. Red Cell Slip at Wall in vitro. Science 1967,15,15541556. Bird, R. B.; Stewart, W. E.; Lightfoot, E. M. Transport Phenomena; John Wiley and Sons Inc.: New York, NY,1960. Caro, C . G.; Fitz-Gerald, J. M.; Schroter, R. C. Atheroma and Arterial Wall Shear: Observation, Correlation, and Proposal of a Shear-dependent Mechanism for Atherogenesis. Proc. R . SOC.London, B 1971,117, 187-205. Dawson, T. H. Engineering Design of the Cardiovascular System of Mammals, 1st ed.; Prentice-Hall: Englewood Cliffs, NJ, 1991. DeBakey, M. E.; Lawrie, G. M.; Glaeser, D. H. Patterns of Atherosclerosis and their Significance. Ann. Surg. 1985,201 (2), 115-131. Denn, Morton M. Process Fluid Mechanics, 1st ed.; Prentice-Hall, Inc.: Englewood Cliffs, NJ, 1980.
Fahraeus, R.; Lindqvist, T. The Viscosity of Blood in Narrow Capillary Tubes. Am. J . Physiol. 1931, 6, 562-568. Fairbairn, J. F.; Juergens, J. L.; Spittell, J. A. Peripheral Vascular Diseases, 4th ed.; W. B. Saunders Company: Philadelphia, PA, 1972. Fry, D. L. Acute Vascular Endothelial Changes Associated with Increased Blood Velocity Gradients. Circ. Res. 1968,22, 165197. Glagov, S.; Newman, W. P.; Schaffer, S. A. Pathobiology of the Human Atherosclerotic Plaque; Springer-Verlag: New York, NY, 1990. Gray, W. G. A Derivation of the Equations for Multi-phase Transport. Chem. Eng. Sei. 1976,30,229-233. Horner, M. Chemical Engineering Honors Thesis, Florida State University, Tallahassee, FL, 1594. Hurst, J. W. The Hearth, 7th ed.; McGrawHill: New York, NY, 1990. Lee, J.; Fung, Y. Flow in Locally Constricted Tubes at Low Reynolds Numbers. J . Appl. Math. 1970,37 (l), 9-16. Liepsch, D. W. Blood Flow in Large Arteries Applications to Atherogenesis and Clinical Medicine; 1st ed.; Karger: New York, NY,1990; p 97. Lightfoot, E. N. Transport Phenomena and Living Systems, 1st ed.; John Wiley & Sons, Inc.: New York, NY,1974; pp 25-26. Merrill, E. W.; Benis, A. M.; Gilliland, E. R.; Sherwood, T. K.; Salzman, E. W. Pressure flow Relations of Human Blood in Hollow Fibers at Low Flow Rates. J . Appl. Physiol. 1965,20, 954-967. Misra, J. C.; Kar, B. K. Momentum Integral Method for Studying Flow Characteristics of Blood Through a Stenosed Vessel. Biorheology 1989,26, 23-25. Nazemi, M.; Kleinstreuer, C. Analysis of Particle Trajectories in Aortic Artery Bifurcations with Stenosis. J . Biomech. Eng. 1989,111,311-315. Nerem, R. M. Vascular Fluid Mechanics, the Arterial Wall, and Atherosclerosis. J . Biomech. Eng. 1992, 114, 274. Paine, M. A.; Carbonell, R. G.; Whitaker, S. Dispersion in Pulsed Systems-I. Heterogeneous Reaction and Reversible Adsorption in Capillary Tubes. Chem. Eng. Sci. 1983,38, 1781. Pedley, T. J. The Fluid Mechanics of Large Blood Vessels, 1st ed.; Cambridge University Press: New York, NY,1980. Perry, R. H.; Green, D. Perry’s Chemical Engineers’ Handbook, 6th ed.; McGraw-Hill Book Company: New York, NY,1984. Schneiderman, G.; Ellis, C. G.; Goldstick, T. K. Mass Transport to Walls of Arteries: Variation with Reynolds Number and Blood Flow Separation. J. Biomech. 1979, 12, 869-877. Slattery, J. C. Momentum, Energy, and Mass Transfer in Continua; Krieger: Melbourne, FL, 1981. Stryer, R. Biochemistry, 3rd ed.; W. H. Freeman and Company: New York, NY,1988; p 561. Taylor, G. I. Dispersion of Soluble Matter in Solvent Flowing Slowly Through a Tube. Proc. R . SOC.London, A 1963, 219, 186-203. Wayland, H. Rheology and the Microcirculation. Gastroenterology 1967,52, 344-346. Whitaker, S. The Spatial Averaging Theorem Revisited. Chem. Eng. Sei. 1985, 40, 1387. Young, D. F.; Tsai, F. Y. Flow Characteristics in Models ofArteria1 Stenoses. J . Biomech. 1973, 6, 395-410. Zarins, C. K.; Giddens, D. P.; Bharadvaj, B. K.; Sottiurai, V. S.; Mabon, R. F.; Glagov, S. Carotid Bifurcation Atherosclerosis: Quantitative Correlation of Plaque Localization with Flow Velocity Profiles and Wall Shear Stresses. Circ. Res. 1983,53 (4), 502-514. Received for review January 5, 1995 Revised manuscript received July 10, 1995 Accepted July 25, 1995@ IE950026C Abstract published in Advance ACS Abstracts, September 1, 1995. @
