Article Outline

Article Information

PubMed

Related Articles

  • …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 6, 1858-1877, 15 March 2007

doi:10.1529/biophysj.106.095042

Biophysical Theory and Modeling

Mesoscale Simulation of Blood Flow in Small Vessels

Prosenjit BagchiGo To Corresponding Author,  

Department of Mechanical and Aerospace Engineering, Rutgers University, The State University of New Jersey, Piscataway, New Jersey

Address reprint requests to Prosenjit Bagchi.

Abstract

Computational modeling of blood flow in microvessels with internal diameter 20–500μm is a major challenge. It is because blood in such vessels behaves as a multiphase suspension of deformable particles. A continuum model of blood is not adequate if the motion of individual red blood cells in the suspension is of interest. At the same time, multiple cells, often a few thousands in number, must also be considered to account for cell-cell hydrodynamic interaction. Moreover, the red blood cells (RBCs) are highly deformable. Deformation of the cells must also be considered in the model, as it is a major determinant of many physiologically significant phenomena, such as formation of a cell-free layer, and the Fahraeus-Lindqvist effect. In this article, we present two-dimensional computational simulation of blood flow in vessels of size 20–300μm at discharge hematocrit of 10–60%, taking into consideration the particulate nature of blood and cell deformation. The numerical model is based on the immersed boundary method, and the red blood cells are modeled as liquid capsules. A large RBC population comprising of as many as 2500 cells are simulated. Migration of the cells normal to the wall of the vessel and the formation of the cell-free layer are studied. Results on the trajectory and velocity traces of the RBCs, and their fluctuations are presented. Also presented are the results on the plug-flow velocity profile of blood, the apparent viscosity, and the Fahraeus-Lindqvist effect. The numerical results also allow us to investigate the variation of apparent blood viscosity along the cross-section of a vessel. The computational results are compared with the experimental results. To the best of our knowledge, this article presents the first simulation to simultaneously consider a large ensemble of red blood cells and the cell deformation.

Introduction

Blood is a multiphase fluid that is primarily made of red blood cells (RBCs), white blood cells, and platelets suspended in plasma. Under normal, healthy conditions, a freely suspended RBC is a biconcave discoid with 8-μm diameter and 2-μm thickness. RBCs constitute ∼40–45% of the total blood volume. Being highly deformable particles, RBCs can easily squeeze through the smallest capillaries having internal diameter less than their characteristic size. The particulate nature of the blood and the deformability of the RBCs determine the overall rheological behavior of the blood.

In vitro studies on blood flow through narrow tubes have revealed complex rheological behavior of blood 1,2,3,4,5,6. In large vessels with internal diameter >500μm, blood behaves as a Newtonian fluid with a constant viscosity. In vessels <500μm, blood behaves as a non-Newtonian fluid. In such vessels, the viscosity of blood depends on the vessel diameter. This behavior is known as the Fahraeus-Lindqvist effect 7. It is characterized by a decrease in the apparent blood viscosity as the vessel diameter decreases below 500μm. The minimum apparent viscosity is reached when the tube diameter is ∼8μm. Upon further decrease in tube diameter, the apparent viscosity increases very rapidly. The physical reason behind the Fahraeus-Lindqvist effect is the formation of a cell-free layer near the wall of the tube 8,9. The layer is devoid of RBCs and has a reduced local viscosity. The core of the tube, on the contrary, is rich with RBCs and has a higher local viscosity 10,11,12,13. The extent of the cell-free layer, which depends on the vessel size and hematocrit, is a major factor that determines the apparent viscosity of blood.

The formation of the cell free-layer is due to a migration of the red blood cells lateral to the mainstream flow and away from the wall of the vessel. The lateral migration arises due to the deformation of the red blood cells 14. As per the theory of viscous fluid mechanics (Stokes flow), a perfectly rigid particle does not migrate away from the wall, but a deformable particle does 15. The rate of migration depends on how easily the particle deforms; an easily deformable particle in a parabolic flow migrates faster toward the center of the vessel than a less deformable particle. The rate of migration also depends on the instantaneous radial location of the particle in the vessel. In a dilute suspension, individual RBCs continuously migrate toward the center of the vessel. In a dense suspension, hydrodynamic interaction between adjacent cells also affects their motion. The cell-free layer is formed under a balance between the deformation-induced lateral migration and the dispersion due to the cell-cell interaction 14.

Significant progress has been made in understanding the mechanical behavior of the red blood cells that have provided a basis for the study of cell deformation. Computational modeling and simulation of blood flow in microvessels have focused mainly on two issues: axisymmetric motion of a single cell in capillary tubes with diameter <10μm 17,18,19, and the deformation of dilutely suspended cells in simple (linear) shear flows 20. Nonaxisymmetric motion of RBCs in cylindrical capillaries has also been addressed 21. Development of more realistic simulation of multiple red blood cells flowing through vessels of diameter 10–500μm has remained a major challenge. It is because blood in such vessels behaves as a multiphase suspension of deformable particles. While a continuum description of blood is sufficient for average flow profile, it is not so if the motion of individual cells and their interaction are concerned. The size of the RBCs is comparable to the size of the vessels, and hence each cell must be taken into consideration in the modeling. At the same time, multiple cells, often of the order of a few thousands in number, must be considered to realistically simulate a microvessel. Any computational model aimed at understanding blood flow in microvessels must take into consideration the deformability of the individual red blood cell, and an ensemble of a large population of cells.

A significant step in the simulations of flowing multiple RBCs has been achieved by Sun and Munn 22, who used a Lattice-Boltzmann simulation to address blood flow in 20–40μm two-dimensional channels and at hematocrit of 10–30%. The simulations were two-dimensional, but the results showed that apparent viscosity of the suspension increases with increasing volume fraction of the particles, which is in agreement with the earlier observations with suspensions of RBC, and other rigid particles. However, the RBCs in their simulation were modeled as rigid disks, rather than deformable particles. Further, the range of channel size is not large enough to show the nonlinear behavior of the apparent viscosity with varying channel size. It should also be noted that at very high volume fractions, which were not considered by Sun and Munn, rigid particles can stop the flow. This is not the case with RBCs due to their flexibility.

In this article, we present computational simulation of the motion of red blood cells flowing through two-dimensional channels of size 20–300μm. Similar to the work of Sun and Munn 22, we consider two-dimensional simulations. However, the deformability of the cells is included in our model. Moreover, a large cell population comprising of as high as 2500 red blood cells are simulated. To the best of our knowledge, this article presents the first simulation to consider such a large ensemble of deformable cells, though in two dimensions. As we will see later, many characteristics of the RBC motion, formation of the cell-free layer, and the Fahraeus-Lindqvist effect are quite accurately predicted by our two-dimensional simulations.

The structure of the article is as follows. In the next section, we describe the simulation technique, followed by presentation of the results. First, we describe the motion of an isolated red blood cell flowing through a vessel, and address its dynamic behavior, such as lateral migration, tank-treading, and flipping motion. We then consider simulations of multiple cells in the range of vessel size 20–300μm and discharge hematocrit 10–60%. Number of RBCs considered in a typical simulation ranges from 5 to 2500. Results on the trajectory and velocity of individual cells and their fluctuation statistics are presented. The inclusion of deformability allows us to study the formation of the cell-free layer. Comparisons are made with the experimental results 1,2,3,4. We then present the apparent viscosity of blood for varying hematocrit and vessel diameter, and discuss the Fahraeus-Lindqvist effect. The numerical results also allow us to investigate variation of “local” apparent viscosity across the cross-section of the vessel.


Simulation technique

Flow configuration

The flow configuration is described in Fig. 1. The motion of an ensemble of red blood cells through a two-dimensional rectangular channel is considered. The undeformed resting shape of a red blood cell is taken to be a biconcave disk. In the figure, the flow is along the x direction and from left to right. The flow is driven by a constant pressure gradient. No-slip conditions are imposed at the walls of the channel. In absence of the cells, the velocity profile of the pure plasma is parabolic and is given by the Poiseuille law. The height of the channel is denoted by H, which is equivalent to the tube diameter in case of a three-dimensional flow. The computation domain is a square segment of sides of length H. Periodic conditions are imposed at the inflow (left) and outflow (right) of the domain. The cells that leave the domain through the outflow boundary are brought back into the domain through the inflow boundary.

Display large version of this figure
Figure 1
Schematic of the computational domain, and the Eulerian and Lagrangian grids.

Immersed boundary method

The simulation technique considered here is the immersed boundary method developed by Peskin 23, and later extended by Tryggvason and co-workers 24,25 as the front-tracking method for deformable interface. The method has been applied to the simulations of suspension of liquid drops and bubbles, and deformation of a red blood cell ghost in a shear flow 20. The method is particularly suitable for this study, as the red blood cells are modeled as liquid “capsules.” The structure of a red blood cell consists of hemoglobin solution surrounded by a lipid bilayer membrane. The liquidlike nature of hemoglobin, and the elastic nature of the membrane give rise to the deformability of the cell. On a mesoscopic scale, the detailed molecular structure of the lipid bilayer and the underlying two-dimensional cytoskeleton network can be neglected. Then, the individual RBC can be modeled as a liquid capsule, that is, a viscous liquid drop surrounded by a thin elastic membrane. The viscosity of the liquid interior (i.e., hemoglobin) of the capsule is five times higher than that of the exterior liquid (i.e., plasma).

In the present model, blood plasma and RBC hemoglobin are assumed to behave as Newtonian fluids. The Newtonian nature of these fluids is well established 26. The non-Newtonian behavior of the whole blood primarily arises due to deformability of individual RBC. The motion of the liquids, plasma, and hemoglobin, is governed by the continuity and Navier-Stokes equations as

(1)
(2)
where u(x, t) is the fluid velocity, ρ is the density, p pressure, and F is a sourcelike term that arises due to the elastic force generated in the cell membrane. The viscous stress tensor τ is given by
(3)
where μ(x) is the viscosity of either the plasma or the hemoglobin solution. For any point within the RBC, μ=μr, and for any point outside the RBC, μ=μp, where μr is the viscosity of the RBC hemoglobin, and μp is the viscosity of the plasma.

The main idea of the front-tracking method is to use a single set of equations for both fluids, plasma and hemoglobin, as in Eq. (2). The Navier-Stokes equations for the fluid flow are solved on a fixed Eulerian grid, and the cell-plasma interface is tracked in a Lagrangian manner by a set of moving grid, used to discretize the cell membrane as shown in Fig. 1. As the cells deform during their motion, the cell membranes are stretched, and elastic forces are generated on the membrane. The deformation of individual cell alters the surrounding flow. The elastic forces at the cell membranes are coupled to the bulk fluid motion via the source term F in Eq. (2) as

(4)
where f(x′, t) is the elastic force generated in the membrane. Here x is the location of an arbitrary point in the flow domain, x′ is any point on the membrane, and δ is the Delta function, which vanishes everywhere except at the membrane. Models for computation of f(x′, t) are described later.

The Navier-Stokes equations are first solved to obtain the fluid velocity and pressure. Then the cells are advected in a Lagrangian manner. The velocity of the cell membrane is obtained by interpolating the velocity of the fluid as

(5)
where S denotes the entire flow domain. The membrane is then advected by
(6)

The δ-functions used in Eqs. (4) and (5) are constructed by multiplying one-dimensional delta functions, such as δ(x-x′)=δ(x-x′)δ(y-y′), in two dimensions. For numerical implementation, a smooth representation of the δ-function is introduced as 24

(7)
where Δ is the Eulerian grid size. The above representation approaches the actual delta function as the grid size approaches zero. The discrete δ-function is so constructed that properties, such as viscosity, vary smoothly over four Eulerian grid points surrounding the RBC membrane. In discrete form, the integrals in Eqs. (4) and (5) can be written as
(8)
(9)
where i and j represent the Lagrangian and Eulerian grid points, respectively. The above representation needs to be modified if the interface is located close to a solid boundary. In our simulation, we ensure that the minimum distance between the solid boundary and the interface is at least four Eulerian grid size.

As the cells move to new positions, the viscosity μ(x, t) needs to be updated. Following Tryggvason and Unverdi and Tryggvason 24,25, this is done by first defining an indicator function I(x) such that

(10)
To find the indicator function, we use a Poisson solver as
(11)
where
(12)
and where D is the discrete δ-function as given in Eq. (7), n is the unit normal to the cell surface, and Δs is a discrete line segment on the cell surface. In this way, the direct solution of the advection equation for viscosity, and hence the smearing of viscosity profile across the interface are avoided. Here we avoid the details of the method, and refer to articles by Tryggvason and co-workers 24,25.


RBC model

As mentioned before, red blood cells in this study are modeled as liquid “capsules,” that is, viscous liquid drops surrounded by elastic membranes. The viscosity of the liquid interior (i.e., hemoglobin) of the capsule is five times higher than that of the exterior liquid (i.e., plasma). This difference is taken into account in the immersed boundary method described above. The deformation of the cells under a dynamic fluid motion generates an elastic force f(x′, t) in the cell membrane. Computation of this force requires a constitutive law for the material of the membrane. Here we assume that the membrane follows the neo-Hookean law. Note that the neo-Hookean law does not strictly represent the behavior of a red blood cell membrane. An RBC membrane is strongly resistant to area dilatation. On the contrary, the neo-Hookean model employed here does allow area dilatation. Membrane models that restrict area dilatation and hence are more accurate for the RBC have been developed 27,28. The present methodology does allow incorporation of these models. The neo-Hookean model is chosen because of its simplicity. An accurate modeling of cell deformation is not the goal of this article. For the present purpose, a cell model that takes into account deformability is sufficient. Indeed, it is shown later that the neo-Hookean model can effectively capture some general characteristics of the RBC motion in a shear flow, such as the tank-treading motion and the lateral migration.

For a two-dimensional neo-Hookean membrane of a three-dimensional cell, the strain energy function is given by 20

(13)
where Es is the shear modulus of elasticity of the membrane, h is the thickness, and ϵ1 and ϵ2 are the principal stretch ratios. The tensions T1 and T2 in the principal directions are then given by 29
(14)
For the two-dimensional simulations considered here, the membrane is a closed curve. A two-dimensional cell is then equivalent to an actual three-dimensional cell subject to a stretching in one direction only. That is, T1 ≠ 0, T2=0, where 1 indicates the in-plane direction along the membrane, and 2 indicates the out-of-plane direction normal to Fig. 1. The deformation ϵ2 in the out-of-plane direction is not zero. But since T2=0, we can express ϵ2 in terms of ϵ1. Then, for a two-dimensional cell, we have
(15)
where T=T1 and ϵ=ϵ1. For a discretized cell, T is the tension acting along a line segment connecting two adjacent Lagrangian grid points on the membrane, and ϵ is the stretch ratio (undeformed length by deformed length) of the line segment. At any Lagrangian grid point on the membrane, two line segments meet. The membrane elastic force f is then the resultant vector of the tensions in the two adjacent segments,
(16)
where i and j denotes two adjacent line segments, and ei and ej are the unit tangent vectors along them.

The RBC membrane also has a bending resistance. To include the bending resistance in our simulation, we follow the approach by Pozrikidis 30. Equation (16) is then modified as

(17)
where q=dm/dl is the transverse shear tension, m=EB(κ-κr) is the bending moment, EB is the bending modulus, κ is the local curvature, κr is the reference curvature in resting configuration, l is the arc length along the membrane, and n is the unit normal vector at a Lagrangian marker point on the cell surface.


Dimensionless parameters, finite difference scheme, and resolution

In absence of the cells, the maximum velocity of the parabolic flow at the channel center is denoted by Ucl. The governing equations are made dimensionless using the channel height H as the characteristic length scale, Ucl as the velocity scale, and H/Ucl as the timescale. In dimensionless form, the shear modulus of elasticity of the RBC membrane is given by E*=μpUcl/Esh, which is the ratio of the viscous force to the elastic force of the capsule membrane. The dimensionless bending stiffness is expressed as EB*=EB/(a2Es), where a is a characteristic dimension of the cell. The Reynolds number of individual RBC, defined as Re=ρUcla/μp, is much less than unity.

In the present model, the dimensionless parameters E* and EB* determine the deformability of the liquid capsule, and hence of the red blood cell. For a normal, healthy red blood cell, Es=0.006 dyn/cm, and EB ≈ 1.8×10−12 dyn-cm 6. Under diseased conditions, e.g., in sickle cell anemia, the cells lose their deformability. The loss of deformability can be expressed in our model in terms of higher-than-normal values of Es and EB. Accordingly, the values of E* and EB* change as the cell loses its deformability. The value of E* also depends on the flow velocity Ucl. Note that in presence of the RBCs, the maximum centerline velocity is significantly reduced below Ucl if the pressure gradient is kept constant (shown later in Plug-Flow Profile). The average flow velocity in presence of the RBCs is denoted by Um. In our simulations, Um ranges from ∼3 mm/s to 15 mm/s (Table 2). Accordingly, the values of E* and EB* used in our simulations ranges from ∼0.02–1.0, and 0.0005–0.002, respectively. The higher values of E* and lower values of EB* typically represent a normal, deformable RBC, whereas lower values of E* and higher values of EB* represent a less deformable RBC. The higher values of E* also represent higher flow velocities, and vice versa. The range of Um considered here matches with that in the in vivo experiments 1,2. A pseudo-shear rate can also be defined as Um/H, which varies between ∼30 and 400s−1, similar to the range in Bishop et al. 1,2. In terms of physical time, the simulations represent 0.1s of flow, on the average.

The governing equations are discretized spatially using a finite difference scheme, and temporally using a two-step time-split scheme. In this method, the momentum equation is split into an advection-diffusion equation and a Poisson equation for the pressure. The body-force term is retained in the advection-diffusion equation. The nonlinear term in this equation is treated explicitly using a second-order Adams-Bashforth scheme. To avoid a stability problem, we treat the viscous terms implicitly using an alternating-direction implicit scheme. In the method, three one-dimensional implicit equations are obtained, which are solved directly by a tridiagonal matrix solver. The velocity is not divergence-free at the end of the advection-diffusion step. The Poisson equation is then solved to obtain pressure at the next time level. Using the new pressure, the velocity field is corrected to make it divergence-free. To reduce expensive computation, the Poisson equation is Fourier-transformed in the periodic direction yielding a set of one-dimensional decoupled PDEs, which is directly inverted to obtain pressure. Details of the time-step scheme are given in Bagchi and Balachandar 31.

The accuracy of the simulations depends on the resolution of the Eulerian and Lagrangian grids. A detail study of the resolution and validation of the computational model are given elsewhere and not repeated here 32. In deciding the resolution, we make sure that there is a sufficient number of Eulerian points within each cell area, and in the region between two adjacent cells. Typically, for a circular cell, ∼25 Eulerian points per diameter are found to be sufficient 33.

The resolution used in the present simulations is given in Table 1. It varies from 129×128 Eulerian grids for a 20-μm channel to 2049×2048 grids for a 300-μm channel. The Lagrangian resolution varies from 128 to 512 marker points per cell. The Lagrangian resolution is increased as the cell deformability increases to ensure that strong curvatures in the cell shape are well resolved. The requirement for a high Eulerian resolution renders some of our computations very expensive, even in two dimensions. Efficient algorithm based on fast Fourier transform, and OpenMP parallelization, have been implemented to speed up the computation. The simulation for a 80-μm channel with 501×500 Eulerian resolution takes ∼50 CPU hours on 1.6GHz IBM p690 processors for 50,000 timesteps.

Table 1 Channel size, discharge, and tube hematocrits, and Eulerian grid resolutions used in the present simulation
Channel size (H, μm)Hd %Ht %Number of RBCEulerian resolution
2020123129×128
2030205129×128
2045337129×128
20604810129×128
40106.45249×248
402013.511249×248
40453528249×248
40605040249×248
80107.632501×500
802015.766501×500
804538160501×500
806054227501×500
15020182511025×1024
15045436001025×1024
15060578001025×1024
300454425002049×2048

A typical dimensionless time-step size used is ∼0.001. In the present method, the immersed boundaries are advected explicitly. In many cases, however, the explicit treatment of the immersed boundary results in more restrictive stability conditions than the viscous terms of the Navier-Stokes equations. For the present simulations, this is not the case due to the specific constitutive law used.



Results

Motion of an isolated RBC

First we describe the motion of a single, isolated RBC in a parabolic flow in a rectangular channel of H=40μm. The results are shown in Figure 2 and Figure 3. The initial resting shape of the cell is biconcave. At time t=0, the cell is located close to the wall of the channel. As the flow starts, the cell deforms and moves longitudinally along the flow direction as well as laterally normal to the flow. Three cases are simulated to study the role of cell deformability: case a, a normal RBC at E*=0.2 and EB*=0.0005; case b, a less deformable RBC with E*=0.02 and EB*=0.002; and case c, a RBC with reduced membrane resistance (E*=1.0, EB*=0.0005).

Display large version of this figure
Figure 2
Motion of an isolated RBC in parabolic flow through a 40-μm channel. The flow is from left to right and the center of the channel is shown by (- · - · - · -). (a) Normal, deformable RBC; (b) less deformable RBC; and (c) RBC with reduced membrane resistance. A point on the cell surface is shown to illustrate the tank-treading motion.
Display large version of this figure
Figure 3
(a) The x-y trajectory, (b) longitudinal velocity, and (c) lateral velocity of the cells shown in Fig. 2. (··· ···) Normal cell (case a in Fig. 2); (——) less deformable cell (case b); (- - - -) cells with reduced membrane resistance (case c).

As shown in Fig. 2, the RBC in cases a and b undergoes significant deformation during its motion in the parabolic flow through the channel. A normal cell, as in case a, repeatedly attains biconcave shape and elliptic shape. The biconcave shape is attained when the major axis of the cell is aligned nearly normal to the flow direction. The elliptic shape is attained when the major axis is between 0° and 45° with the flow direction. At reduced membrane resistance, as in case c, the biconcave shape is completely lost, and the cell attains a nearly elliptic shape. For the less deformable cell, as in case b, no significant deformation is observed, and the cell maintains the biconcave shape throughout its motion.

For an RBC placed in a shear flow, two modes of motion have been observed by previous researchers, both experimentally and computationally: tank-treading motion and tumbling motion 14,18,34,35. In the tank-treading mode, the cell membrane and the interior liquid undergo steady rotary motion while the cell maintains a fixed orientation with the flow. In the tumbling motion, the cell flips like a rigid body. The transition from tank-treading to tumbling motion occurs as the deformability of cell decreases. Our results in Fig. 2 reproduce these earlier observations. For the RBC with reduced membrane resistance, as in case c, only tank-treading motion is observed. In the figure, an arbitrary point on the cell surface is marked to show the tank-treading motion. For the less deformable RBC, as in case b, only tumbling motion is observed. On the contrary, for the normal RBC, as shown in case a, simultaneous tumbling and tank-treading motions are observed. The frequency of tumbling motion increases as the deformability decreases.

For all cases considered, the RBC is observed to migrate laterally away from the wall toward the center of the channel under the action of the parabolic flow. The lateral position of the center of the RBC is shown in Figure 3a. In general, the migration is a very slow process; for the normal RBC, the cell travels only 10μm in the lateral direction while moving nearly 2500μm in the longitudinal direction. The rate of migration depends on the deformability of the cell. Rate of migration is lower for the less deformable RBC as in case b than in cases a and c in Fig. 2. Interestingly, we also observe that the migration rate is higher for the normal cell, as in case a, than that for a cell with reduced membrane resistance, as in case c. Note that the normal RBC in case a performs both deformation and tumbling motion. On the contrary, no tumbling motion is seen in case c. The results suggest a possible coupling between deformation and the tumbling motion, which results in a higher migration rate in case a.

The longitudinal and the lateral velocity components of the RBC are shown in Figure 3bc. Clearly, the lateral component is an order-of-magnitude less than the longitudinal component. The velocity components, as well as the lateral position, show fluctuations for the normal cell and the less deformable cell. These fluctuations arise due to the tumbling motion. The lateral velocity may become periodically negative due to the tumbling motion as in the case of the less deformable cell. Within one cycle of oscillation, the migration velocity becomes maximum when the cell is aligned nearly at 45° with the flow, and minimum at 135°. Oscillation increases with decreasing deformability.

The above results represent the dynamics of a red blood cell in a dilute suspension flowing through a conduit. The results presented here on the lateral migration of the RBC agree with the glass tube experiments 14. Clearly, the present computational model is able to capture the general dynamic behavior of a red blood cell in a parabolic flow, particularly the tank-treading and tumbling motion, and the lateral migration. As mentioned before, the lateral migration leads to the formation of the cell-free layer, which serves as the primary mechanism for the Fahraeus-Lindqvist effect. In a nondilute suspension, the presence of many RBCs affects the motion of individual cell, and the formation of the cell-free layer. It is of interest in the next section to study how the motion of individual cell is affected in presence of neighboring cells.

It should be noted that according to Keller and Skalak 35, an ellipsoid with an internal-to-external viscosity ratio of 5 would tumble rather than rotate even at high shear rate. The results in Figure 2c may appear to be in contrast to their result. In Keller and Skalak 35, the particles are nondeformable. For a deformable cell, the transition from flipping to tank-treading motion also depends on the aspect ratio of the cell, and hence the extensional resistance of the membrane. Note that the neo-Hookean model used here for the cell membrane does allow continuous extension of the membrane with increasing shear rate. Ramanujan and Pozrikidis 36 considered the large-deformation of ellipsoidal and biconcave capsule in shear flow using neo-Hookean model, and showed that at viscosity ratio of 5, an ellipsoidal cell performs only oscillatory motion rather than a flipping motion. Their result indicates that the ellipsoid is eventually likely to achieve a steady orientation. The biconcave discoid however showed a flipping motion. If the initial shape is spherical, the deformed ellipsoidal cell does not show even an oscillatory motion at viscosity ratio 5.

As far as the deformation of a cell is concerned, the viscosity ratio and the elastic resistance of the cell membrane contribute in the same way. That is, an increase in any of these two parameters would cause less deformation. Since high viscosity ratio causes transition from tank-treading to flipping motion, so likely does the higher membrane resistance. Thus allowing the membrane deformability, which is neglected in Keller and Skalak 35, may delay the transition from tank-treading to flipping motion only to a viscosity ratio >5. In fact, if membrane deformability is allowed, as done in our article, the cell would elongate more at a given viscosity ratio. Keller and Skalak 35 mentioned that increasing the elongation promotes a stationary orientation of the cell. Thus, the result in Figure 2c is not completely in opposite to that of Keller and Skalak 35.


Motion of RBC suspension

We now present the results on the simulation of suspension of multiple red blood cells. As mentioned before, the size of the vessel ranges from H=20–300μm, and the discharge hematocrit Hd=10–60%. The number of red blood cells considered in our simulations varies from 5 to 2500. Note that in the simulations, discharge hematocrit is not directly specified. Instead, we specify the tube hematocrit Ht that varies from 6 to 57%. In Table 1, we have listed the discharge hematocrit, corresponding tube hematocrit, and the number of RBCs in various numerical experiments considered here. In the subsequent results, both tube and discharge hematocrits are mentioned. For a given discharge hematocrit, the tube hematocrit is first obtained as 3

(18)
where D≡H is the tube diameter. Once the tube hematocrit is obtained, the number of red blood cells within the computation domain is found by knowing the cell volume (or, area, in two dimensions). The computation domain is a square segment of the channel. The initial shape of the cells is biconcave, and the cells are initially distributed in a random manner throughout the domain. The flow starts at time t=0 under a constant pressure gradient. In the simulations, the velocity and pressure fields of the fluid, and the coordinates, shapes, and velocity of the cells are stored at frequent time intervals. Instantaneous distributions of the red blood cells in the vessels are shown in Figure 4 and Figure 5 and Figure 6 and Figure 7 for a few representative cases. Simulations are performed over sufficient time so that quasi-steady state is reached.

Display large version of this figure
Figure 4
RBC suspension in a 20-μm channel. (a) Normal RBC at Ht=20% (Hd=30%); (b) normal RBC at Ht=48% (Hd=60%); and (c) less deformable RBCs at Ht=20% (Hd=30%). For each case, three time instances are shown. The mean velocities are (a) 7.5 mm/s; (b) 3.5 mm/s; and (c) 4.7 mm/s. The third plots for panels a–c represent ∼500ms after the onset of flow.
Display large version of this figure
Figure 5
Suspension of normal RBCs at Hd=45% in (a) 80-μm channel, and in (b) 150-μm channel. The tube hematocrits are 38% and 43%, respectively. The computation domain contains 160 cells in panel a and 600 cells in panel b. Mean velocities are (a) 5 mm/s, (b) 6.5 mm/s. The figures represent ∼300ms after the onset of flow.
Display large version of this figure
Figure 6
Suspension of normal RBCs in 300-μm channel at Hd=45% (Ht=44%). A total of 2500 cells are simulated in the computation domain as shown above. Mean velocity is 12 mm/s. The figure represents ∼150ms after the onset of flow.
Display large version of this figure
Figure 7
Suspension of less deformable RBCs at Hd=45% in (a) 80-μm channel and (b) 150-μm channel. The tube hematocrits are 38% and 43%, respectively. The computation domain contains 160 cells in panel a and 600 cells in panel b. Mean velocities are (a) 3.6 mm/s and (b) 4.9 mm/s. The figures represent 300ms after the onset of flow.

Shown in Fig. 4 are the results for a 20-μm channel. Three different cases are considered here: suspension of 1), normal RBCs at Ht=20% (Hd=30%); 2), normal RBCs at Ht=48% (Hd=60%); and 3, less deformable RBCs at Ht=20% (Hd=30%). Time evolution of the cell distribution is shown in the figure, and a few cells are marked by numbers. First consider case a, for normal RBCs at Ht=20%. As the flow develops, the cells migrate toward the center of the channel, and the regions near the walls become devoid of RBCs. However, continuous lateral migration is prevented due to the presence of the neighboring cells. A balance between the hydrodynamic interactions among the cells and lateral migration of individual cell is attained, and a cell-free layer near the wall develops. Significant deformation of the RBCs is observed. Nearly all cells lose their biconcave shape. Unlike the case of a single, isolated RBC as considered in the previous section, none of the cells in suspension performs tumbling motion. The repeated emergence of the biconcave and elliptic shapes as observed before are also not seen here. However, the shapes are changing continuously due to the interaction with the neighboring cells. The RBCs near the center assume slipper shapes, whereas those further away from the center assume nearly elliptic shapes. The slipper shapes of the RBCs in vessels of this range of size have been observed in experiments and previous numerical simulations 6,37.

When the tube hematocrit is increased to Ht=48% (Hd=60%), the cells are more evenly distributed across the channel. The cell-free layer is not well developed. The continuously changing shapes of the cells are still observed, although the slipper shape near the center is less common now. For less deformable RBCs at Ht=20%, as in case c, the initial biconcave shape of individual cell is retained through the simulation. The tumbling motion of the cells is evident here. Tumbling motion results in a higher dispersion of the cells. As a result, the cell-free layer near the wall is reduced compared to that for the normal RBCs at the same Hd. The interface between the cell-free and cell-rich regions is also not well defined due to the cell-cell interaction.

The results for normal RBCs in 80- and 150-μm channels are shown in Fig. 5, while those in 300-μm channel are shown in Fig. 6. In these figures, discharge hematocrit is kept constant at Hd=45%. The tube hematocrits are 38%, 43%, and 44%, respectively. As the vessel size increases, the slipper shape of individual RBC is no longer observed. However, deformation of red blood cells is evident in all cases. The figures show that the cells near the wall deform more and lose their biconcave shape, whereas the cells near the center deform less and retain the biconcave shape. This behavior is expected, since the fluid shear rate decreases from the wall toward the center. Tumbling motion of the cells appears to be suppressed. Most cells near the wall are aligned at an angle with the flow direction, whereas the cells near the center are nearly vertical or parallel to the flow.

The results for less deformable cells are shown in Fig. 7 for 80- and 150-μm channels. As expected, less deformable cells retain their biconcave shape, and perform the tumbling motion. The tumbling motion is stronger for the cells located near the wall than those located near the center. This observation is consistent with the hydrodynamic theory of particle motion in shear flow. The rate of tumbling is proportional to the shear rate of the fluid, which decreases from the wall toward the center of the vessel. However, due to higher hematocrit, tumbling motion is not as strong as that of an isolated RBC as was observed in Figure 2b. In general, the direction of the tumbling motion of individual cells in suspension matches with that of an isolated cell. In the lower half of the vessel, the cells tumble in the clockwise direction, whereas in the upper half they tumble counterclockwise, in accordance with the direction of vorticity of the flow. However, in some cases, strong cell-cell interaction is observed to reverse the direction of rotation. The tumbling motion combined with the cell-cell interaction results in random orientation of the cells across the channel.

It should be noted that in experiments with flowing RBC suspension, the biconcave shape is usually not observed. In the present simulations, we have considered both normal and hardened cells. The hardened cells, as considered in Fig. 7, are expected to retain the biconcave shape due to the high bending-resistance values used to model them. The normal cell, after sufficient simulation time, would lose its biconcave shape. This is evident in Figure 4ab. In Fig. 5, most cells near the channel wall lose the biconcave shape due to local high shear rate. However, these cells are tumbling also. So their shape repeats between biconcave and elliptic shapes, as it was observed for an isolated cell in Figure 2a. As for Fig. 6, the simulation is performed for a short time, since this case is computationally expensive. The two-dimensionality of the problem can also affect the cell shape as a real three-dimensional cell would deform more easily and may lose the biconcave shape on a shorter timescale than a two-dimensional cell.


RBC trajectory and velocity traces

In our simulations, the position and velocity of all red blood cells are tracked in time. These data allow us to study the trajectory and instantaneous velocity of individual red blood cell in the suspension. The trajectory of a few cells are shown in Figure 8 and Figure 9, and the velocity traces are shown in Fig. 10. The cells were tracked over a longitudinal distance that ranges from 500 to 2000μm, depending on the specific simulation. As evident in the figures, individual red blood cells exhibit fluctuations in lateral position and velocity. Fluctuations arise due to the tumbling motion of individual cells as well as from interaction with neighboring cells.

Display large version of this figure
Figure 8
Trajectory of normal red blood cells in suspension flowing through a vessel. Panels a–c represent a 40-μm channel at Hd=10, 20, and 60%, respectively. Corresponding tube hematocrits are 6.4, 13.5, and 50%. (d) 150-μm channel at Hd=20% (Ht=18%). For panel d, only half of the channel is shown.
Display large version of this figure
Figure 9
Trajectory of red blood cells in suspension flowing through a vessel. (a) Normal cells in a 150-μm channel at Hd=45% (Ht=43%); (b) normal cells in a 300-μm channel Hd=45% (Ht=44%); (c) less deformable cells in a 40-μm channel at Hd=20% (Ht=13.5%); and (d) fewer deformable cells in a 40-μm channel at Hd=60% (Ht=50%). In panels a and b, only half of the channel is shown.
Display large version of this figure
Figure 10
Velocity traces of red blood cells in suspension flowing through a 40-μm channel. (a) Normal cells at Hd=20% (Ht=13.5%); (b) normal cells at Hd=60% (Ht=50%); and (c) less deformable cells at Hd=20% (Ht=13.5%). Mean velocity of whole blood for various cases can be found from Table 2. The velocity of RBCs is scaled with the centerline velocity.

First consider a 40-μm channel with normal RBCs, for which three different hematocrits are considered in Fig. 8: case a, Hd=10%; case b, Hd=20%; and case c, Hd=60%. The corresponding tube hematocrits are 6.4, 13.5, and 50%. In case a, lateral migration of the red blood cells initially located close to the wall is observed. However, the rate of migration of a cell in the suspension is much lower than that of an isolated RBC possibly due to cell-cell interaction. Moreover, unlike an isolated cell, the cells in suspension do not migrate continuously. Rather, the trajectories show random fluctuations due to the cell-cell interaction. A comparison of the three cases shows that the fluctuations in the RBC trajectory depend on the hematocrit. For case a at Hd=10%, oscillations in the trajectory are similar to those observed previously in Motion of an Isolated RBC for an isolated red blood cell. Such small-amplitude, low-frequency fluctuations are due to the tumbling motion of individual cells. The cell-cell interaction is less in this case of low Hd. As Hd increases to 20%, as in case b, the trajectory becomes more erratic. Large amplitude but less frequent fluctuations are observed which are due to the increased interaction between the cells. Upon further increase of Hd to 60%, small amplitude frequent fluctuations are nearly suppressed. At this higher hematocrit, the cells move in a nearly stacklike manner. The tumbling motion of individual cells is nearly inhibited, and the fluctuations result mostly from the cell-cell interaction.

The effect of increasing channel size is shown in Figure 9d where the 150-μm channel at Hd=20% (Ht=18%) is considered. Note that for this vessel, only half of the cross section is shown. A slow migration of the cells away from the wall is observed. Fluctuations in the lateral position of the cells indicative of the tumbling motion of individual RBC and cell-cell interactions are observed. In Figure 9a we show the effect of increasing hematocrit while channel size is kept constant at 150μm. The fluctuations in the lateral position appear to diminish at Hd=45% (Ht=43%). Finally, for the 300-μm channels as shown in Figure 9b for Hd=45% (Ht=44%), fluctuations in the trajectory are significantly reduced, and the cells appear to move in a stacklike manner.

The results for less deformable RBCs in a 40-μm channel are shown in Figure 9cd, for Hd=20 and 60%, respectively. The tube hematocrits are 13.5 and 50%, respectively. The trajectory at 20% hematocrit now shows more erratic behavior compared to those for the normal RBCs. This behavior can be explained based on the results presented in the previous sections. A less deformable cell performs a strong tumbling motion, which can significantly affect the motion of the neighboring cells, resulting in more chaotic trajectory. Interestingly, some cells are seen to move toward the wall of the vessel, rather than the center. This anomalous behavior is due to the dispersion of the RBCs resulting from a strong cell-cell interaction. At higher hematocrit, as shown in Figure 9d, the tumbling motion is again inhibited, and the cells appear to move in nearly straight lines as in the case of normal RBCs.

The velocity traces of the RBCs in a 40-μm channel are shown in Fig. 10 for three cases: case a, normal RBCs at Hd=20% (Ht=13.5%); case b, normal RBCs at Hd=60% (Ht=50%); and case c, less deformable RBCs at Hd=20% (Ht=13.5%). In the figure, cells with higher velocity are flowing closer to the center. Oscillations in the velocity traces are larger than those observed earlier for isolated red blood cell. The increased oscillation is due to the cell-cell interaction in the suspension. The oscillations significantly increase as the cells lose deformability, implying increased dispersion due to cell-cell interaction.


Statistics on RBC motion

Root mean-square (RMS) of the fluctuations in the lateral position, and the coefficient of variation (CV) of velocity can be obtained for each red blood cell in the simulations. These quantities are defined as

(19)
and
(20)
where y(t) and v(t) are the lateral location and velocity of the center of a red blood cell at any time instant, and are their mean, and T is the time window over which data is collected. Typically, the RMS and CV are computed over a time in which the cells travel a longitudinal distance of 500–2000μm. Averages on cell statistics are often done using >300 instantaneous measurements.

RMS fluctuations in lateral positions of the red blood cells are shown in Fig. 11 for 40- and 80-μm channels, and for Hd=10–45%. The corresponding tube hematocrits are also mentioned in the figure. The numerical results are compared with the in vivo results 2. In Bishop et al. 2, the velocity and positions of red blood cells were measured in venules (45–75μm diameter) of the rat spinotrapezius muscle. The numerical results in the figure are shown using various symbols corresponding to different vessel size and hematocrit. For the results of Bishop et al. 2, only the range of their data is indicated. The range of RMS lateral positions reported in Bishop et al. 2 is ∼1.5–2.5μm. In comparison, our numerical results yield a range of 0.5–3μm. We also observe that the RMS values show a weak dependence on the vessel size. The RMS of fluctuations in the 80-μm channel is higher than that in the 40-μm vessel. While no universal pattern in the RMS fluctuation is observed, the results for H=40μm and Hd=10 and 20% show that the RMS is higher near the center of the channel, and lower near the wall. At 45% hematocrit, on the contrary, the RMS near the center is lower than that near the wall.

Display large version of this figure
Figure 11
RMS fluctuations in lateral position of red blood cells in suspension. Symbols are present numerical simulations. Range of data obtained from the experiments 2 is also shown.

The CV of velocity is plotted in Fig. 12. Here also we compare the simulation results with the in vivo data 2. The simulation results show the similar pattern as observed in the in vivo data: the value of the CV is higher near the wall of the vessel, and lower near the center. The in vivo data for Hd=45% is in the range 12–22%, while the simulation data is in the range 8–16%, and 4–8% for Hd=20% and 60%, respectively. The CV decreases with increasing hematocrit due to the reduced tumbling motion and cell-cell interaction in an increasingly close-packed arrangement. In Fig. 12, we also show the CV for the less deformable cells, which yield a higher range of 18–24%. This result is expected based on our earlier observation on the velocity traces for the less deformable RBCs.

Display large version of this figure
Figure 12
Coefficient of variation of red blood cell velocity in a 40-μm channel. Symbols are the present numerical simulations. Thick lines are the best fit through the numerical data. Thin line is the best fit from experimental results 2.

Plug-flow profile

Due to the random motion of the red blood cells, the velocity of the bulk fluid is constantly changing in time. The velocity data over the entire computational domain is stored at frequent intervals during the simulation. They are post-processed to obtain the average velocity profile of blood. The mean velocity profile is obtained by averaging nearly 300 instantaneous measurements, and over all grid points along the x direction.

The time-averaged velocity profiles are shown in Fig. 13 for 20–300-μm channels. The mean velocity (averaged over the cross-section) in dimensional form for various cases is given in Table 2. Also shown is the parabolic profile of the Poiseuille flow, which occurs in pure plasma in absence of the cells for the same pressure gradient. The effect of vessel size, hematocrit, and RBC deformability on the velocity profile is studied here. Consider first the 20-μm channel in Figure 13a, for which three different discharge hematocrits (20, 45, and 60%) are shown. The values of the tube hematocrits can be found from Table 1. In presence of the RBCs, the plug-flow profile can be seen which is characterized by a nearly constant velocity near the center of the channel. As discharge hematocrit increases, the plug-flow profile becomes more prominent, and extends toward the wall. The maximum centerline velocity rapidly decreases with increasing hematocrit. Also shown is the velocity profile of the blood with less deformable RBCs. Significant reduction is observed in the centerline velocity compared to that with the normal RBCs at the same Hd.

Display large version of this figure
Figure 13
Average velocity profile of blood. Dotted line is the parabolic flow. Thick lines represent suspension of normal, deformable RBCs, and thin lines represent suspension of less deformable RBCs. (– – – –) Hd=20%; ( – · – · –) Hd=45%; and (– · · – · · –) Hd=60%. In plot d, the solid line represents the 300-μm channel at Hd=45%. The corresponding tube hematocrit can be found from Table 1. Mean velocities are listed in Table 2.
Table 2 Mean blood velocity and shear rate in some representative cases considered in the simulations
Channel size (H, μm)Hd %RBC typeMean velocity (mm/s)Shear rate (s−1)
2030Normal7.5375
2045Normal5.7285
2060Normal3.5175
2030Rigid4.7235
4020Normal13325
4045Normal9.5237
4060Normal7.2180
4020Rigid7.9197
4045Rigid6.5163
4060Rigid6.0150
8020Normal7.594
8045Normal5.063
8060Normal3.240
8045Rigid3.645
15020Normal10.268
15045Normal6.543
15045Rigid4.933
30045Normal1240

The results for 40- and 80-μm channels are shown in Figure 13bc, respectively. Again, the effect of hematocrit and cell deformability is shown. Overall similar behavior of the mean velocity is seen here as in the 20-μm channel. The centerline velocity decreases with increasing Hd and decreasing cell deformability. In the 40-μm channel, the centerline velocity in presence of less deformable cells at Hd=60% is only 20% of Ucl. However, unlike the 20-μm channel, the plug-flow profile is not very prominent for the larger channels. For the 40-μm channel, the profiles are blunt with reduced centerline velocity for Hd=20%. A plug-flow-like profile is observed only at higher Hd. Same is the case for 80-μm channel; at low Hd, the velocity profile is nearly parabolic, and at higher Hd, it is plug-flow type.

The results for 150- and 300-μm channels are shown in Figure 13d. The plug-flow profile is clearly absent in these larger vessels. The velocity profiles appear to be nearly parabolic, but with significantly reduced centerline velocity. The numerical observation here is in agreement with the in vivo observation 1,2 and analytical prediction 13. In Bishop et al. 1,2, a fully plug-flow profile was not seen for normal blood in 45–75μm venules at Hd=40–50%. However, the profiles were blunt with reduced centerline velocity as observed here. The in vivo data 1,2 were fit with the equation u/Umax=1-(r/R)K. Here the exponent K=2 gives a parabola. When K>3, a clear plug-flow profile is obtained. In the experiments 1,2 for normal blood, K lies between 2.1 and 2.2 indicating a blunt profile with reduced centerline velocity. The present numerical data were fit with u/Umax=1-(r/R)K. It is observed that the numerical data also produce the value of K in the same range as in Bishop et al. 1,2. Sharan and Popel 13 showed that the bluntness of the velocity profile decreases with increasing vessel size and decreasing Hd, in agreement with our numerical results.


Cell-free layer

As discussed earlier, the major factor that contributes to the Fahraeus-Lindqvist effect is the formation of a cell-free layer near the wall of a vessel. Since deformability of the cells is taken into consideration in our computational model, the formation of the cell-free layer can be directly studied. Knowledge about the thickness of the cell-free layer is important, since in many two-phase models of blood flow it is taken as an empirical constant. Note that in our simulations, the red blood cells are initially distributed in a random manner throughout the vessel. As the flow starts, individual cells migrate away from the wall of the channel due to the effect of hydrodynamic shear. However, interaction from the neighboring cells also affects their motion. The cells are repeatedly dispersed toward the wall by such interaction. Eventually a quasi-steady state is reached, and the cell-free layer is formed under a balance of the shear-induced migration and cell-cell interaction.

To obtain the cell-free layer, we first calculate the cell number density distribution across the channel. Following Durlofsky and Brady 38 and Zhou and Pozrikidis 39, we divide the channel into several horizontal zones, count the number of cells in each zone, and normalize by the total number of cells. The instantaneous results are then averaged over time. The interface between the cell-free and cell-rich layers is identified as when the number density reaches zero closest to the wall. The thickness of the cell-free layer as obtained from our simulations is presented in Fig. 14. Here the dimensionless thickness δ/(H/2), where δ is the dimensional thickness, is plotted. The effect of vessel size, hematocrit, and cell deformability on the cell-free layer thickness is presented. The in vitro data from Bugliarello and Sevilla 4 and the analytical prediction by Sharan and Popel 13 are also presented, and they show good agreement with the numerical data. The dimensionless thickness of the cell-free layer decreases with increasing vessel size, and increasing hematocrit. For a 20-μm channel with Hd=20%, the cell-free layer covers almost half of the channel. For a 150-μm channel at Hd=10%, the layer covers only 10% of the cross section. A significant decrease in the cell-free layer is observed for blood with less deformable RBCs. At 20% hematocrit, the layer occupies only ∼10% of the cross section in a 40-μm vessel.

Display large version of this figure
Figure 14
Dimensionless cell-free layer δ/(H/2). Lines are analytical modeling 13, and open symbols are in vitro data 4. Solid symbols are the results from the present numerical simulation.

Apparent viscosity of blood

Since in a microvessel, blood does not behave as a Newtonian fluid, the viscosity of the whole blood is expressed in terms of an apparent viscosity, which is defined as

(21)
where D is the tube diameter, Q is the volumetric flow rate, and ΔP/L is the pressure gradient. The apparent viscosity depends on the hematocrit, and vessel size. The Fahraeus-Lindqvist effect refers to the decrease in μapp as the vessel size decreases from 500 to 10μm. In dimensionless form, a relative apparent viscosity is expressed as
(22)
where Qp is the flow rate of the parabolic (Poiseuille) flow in absence of the cells.

The relative apparent viscosity computed from the present simulations is shown in Fig. 15 for the diameter range 20–300μm, and hematocrit range Hd=20–60%. The results from our simulations are compared with the experimental results given in Pries et al. 3. Based on a comprehensive database for viscosity of blood in narrow glass tubes, Pries et al. 3 gave an empirical expression that takes into account both the effect of vessel diameter and hematocrit. Relative viscosity obtained using their expression is shown in the figure as three solid lines representing Hd=20, 45, and 60%. The numerical results are shown using symbols, which agree very well with the experimental fit of Pries et al. 3. The nonlinear decrease in μrel as the vessel size drops from 300 to 20μm is correctly reproduced by our simulations. The increase in μrel with increasing Hd is also correctly predicted by our simulations. It should be mentioned that μrel in our simulations appears to be sensitive, though weakly, to the mean flow velocity Um. In the simulation Um is not specified a priori, rather it is obtained posteriori. Due to the sensitivity of μrel to Um, the actual number of simulations performed is much more than the number of data shown in the figure to obtain a closest match with the empirical fit of Pries et al. 3.

Display large version of this figure
Figure 15
The Fahraeus-Lindqvist effect: the relative viscosity of blood as a function of vessel size and discharge hematocrit. The three solid lines represent the empirical expression given in Pries et al. 3 based on in vitro data. The symbols are the data from the present numerical simulation. The tube hematocrit values for a channel size are given in Table 1.

Local variation of μrel

Flow of blood in microvessels is often described by a two-phase model 13,40,41. In such models, the tube is divided into two regions: a cylindrical core region, and an annular cell-free region. The viscosity of the cell-free layer is usually taken to be equal to that of the plasma. It is noted in Sharan and Popel 13 that the interface between the cell-rich core region and the cell-free layer is not smooth. It is rather rough due to the presence of the red blood cells. Sharan and Popel 13 showed that the apparent viscosity in the cell-free layer is higher than the plasma viscosity. They hypothesized that this was due to the rough interface.

In reality, the interface between the core region and the cell-free layer may not be well defined. The RBCs are continuously dispersed toward the wall due to the hydrodynamic interaction with neighboring cells, as it was observed in the present simulations. The interaction also results in heterogeneous cell distribution. Thus the “effective” viscosity of blood, defined as μ(y)/μp, is expected to vary across the cross-section of the vessel. Damiano 10 used a semiempirical model in which the effective viscosity is assumed to decrease from a maximum value in the core to a minimum value near the wall. The effective viscosity near the wall corresponds to the plasma viscosity. Subsequently, μPIV measurements 11,12 provided experimental verification on such a “local” variation of viscosity across the blood vessel.

The cross-sectional variation of viscosity can be extracted from the present simulations following the approach described in Long et al. 12 and Damiano et al. 11. For a two-dimensional flow of viscous fluid in a channel, with or without cells, the axial momentum equation can be simplified as

(23)
where τxy is the shear stress. Since the pressure gradient is constant, the above equation can be integrated to give
(24)

The constitutive relationship between the local shear stress and rate of strain can then be invoked as

(25)
where and μ(y) are the shear rate and blood viscosity that vary along the cross-section of the channel. The above two relations give
(26)
which is used to obtain μ(y) in our simulations.

Simulation results on μ(y)/μp are shown in Fig. 16. First we note that μ(y) is equal to the plasma viscosity very close to the wall, and it is higher than the plasma viscosity near the center of the vessel. This behavior is in agreement with the assumptions made in the macroscopic two-phase models of blood as mentioned above. However, over the cross section of the vessel, μ(y) shows a strongly nonmonotonic behavior, unlike the monotonic variation assumed in Damiano 10. Most strikingly, μ(y)