# An Anisotropic Constitutive Equation for the Stress Tensor of Blood Based on Mixture Theory

Authors: Mehrdad Massoudi; James F. Antaki
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

Based on ideas proposed by Massoudi and Rajagopal (M-R), we develop a model for blood using the theory of interacting continua, that is, the mixture theory. We first provide a brief review of mixture theory, and then discuss certain issues in constitutive modeling of a two-component mixture. In the present formulation, we ignore the biochemistry of blood and assume that blood is composed of red blood cells (RBCs) suspended in plasma, where the plasma behaves as a linearly viscous fluid and the RBCs are modeled as an anisotropic nonlinear density-gradient-type fluid. We obtain a constitutive relation for blood, based on the simplified constitutive relations derived for plasma and RBCs. A simple shear flow is discussed, and an exact solution is obtained for a very special case; for more general cases, it is necessary to solve the nonlinear coupled equations numerically.


---

## Body

## 1. Introduction

Amongst the  multiphase (more accurately multicomponent) flows
occurring in nature, one can identify blood, mud slides, avalanches, and so
forth, and amongst the many flows found in the chemical industries, one can
name fluidization and gasification, and in agricultural and pharmaceutical
applications, one can name the transport, storage, drying of grains, and so
forth. In most applications, the phases are not of the same material and thus
it is more appropriate and accurate to refer to these cases as “multicomponent”
problems. Historically, two distinct approaches have been used in modeling
multicomponent flows (for simplicity in the remainder of the paper, we assume a
two-component flow). In thefirst case, the amount of the dispersed component is so small that the motion of this
component (measured by a field variable such as concentration) does not greatly
affect the motion of the continuous component (the other component). This is
generally known as the “dilute phase approach,” sometimes also called the Lagrangean approach. This method is used
extensively in applications such as atomization, sprays, and in flows where
bubbles, droplets, and particles are treated as the dispersed rather than a
continuous component. In the second approach, known as the “dense phase approach,” sometimes also called the Eulerian (or the two-fluid) formulation,
the two components interact with each other to such an extent that each
component directly influences the motion and the behavior of the other component.
This method is used extensively in fluidization, gas-solid flows, pneumatic
conveying, suspensions, polymeric solutions, and so forth, (see Kaloni et al.
[1], Massoudi [2] and references therein).The large number of articles published concerning two-component
flows typically employs one of the two continuum theories developed to describe
such situations: mixture theory (the theory of interacting continua) (Rajagopal
and Tao [3]) or Averaging Method(s) (Ishii [4]). Both approaches are based
on the underlying assumption that each component may be mathematically
described as a  continuum. The averaging method directly modifies the classical
transport equations to account for the discontinuities or “jump” conditions at
moving boundaries between the components (cf., Anderson and Jackson [5],
Drew and Segal [6], Gidaspow [7], Jung et al. [8]). The modified
balance equations must then be averaged in either space or time (hence the name
averaging) to arrive at an acceptable local form. In this approach point-wise
equations of motion, valid for a single fluid or a single particle, are
modified to account for the presence of the other component and the
interactions between components. These equations are then averaged over time or
some suitable volume that is large compared with a characteristic dimension
(e.g., particle spacing or the diameter of the particles) but small
compared to the dimensions of the whole system. From the mathematical
manipulation of the averaged quantities, a number of terms, some of unknown
physical origin, arise. These terms are usually interpreted as some form of
interaction between the constituents. Although the two methods seem similar,
the way they approach the formulation of constitutive models is very different.
In fact, as shown in Massoudi [2], many of the interaction models used by
researchers in the averaging community are not frame-indifferent, thus
violating basic principles in physics. Other differences between the two
approaches are explained in Johnson et al. [9, 10], and Massoudi et al. [11, 12].Mixture theory, or the Theory
of Interacting Continua, traces its origins to the work of Fick (1855) (see
Rajagopal [13]) and was first presented within the framework of continuum
mechanics by Truesdell [14]. It is  a means of generalizing the equations and
principles of the mechanics of a single continuum to include any number of
superimposed continua. In an important paper, Rajagopal [13], within the
confines of mixture theory, outlines a procedure to obtain a hierarchy of
approximate models, used in multicomponent situations, including those of Fick
and Darcy. mixture theory is in a sense a homogenization approach in which each
component is regarded as a single continuum and at each instant of time, every
point in space is considered to be occupied by a particle belonging to each
component of the mixture (cf., Truesdell [15]). It provides a means for
studying the motions of bodies made up of several constituents by generalizing
the equations and principles of the mechanics of a single continuum and in
recent years it has been applied to a variety of applications such as fluid-solid particles,
lubrication
with binary-mixtures of bubbly oil, viscoelastic
porous mixtures, swelling  porous
media with microstructure, reacting immiscible mixtures, polymeric solutions,  growth and remodeling of soft tissues,
and ionized fluid mixtures (see Massoudi [16]). More detailed information,
including an account of the historical development, is available in the
articles by Atkin and Craine [17, 18], Bowen [19], Bedford and Drumheller
[20], and in the books by Truesdell [15], Samohyl [21], and Rajagopal and
Tao [3]. mixture theory has also been used in a variety of biomechanics
applications (see, e.g., Ateshian et al. [22], Garikipati et al. [23],
Humphrey and Rajagopal [24], Klisch and Lotz [25], Lemon et al. [26],
Tao et al. [27], Ramtani [28]).It is known that in large vessels (whole) blood behaves as a Navier-Stokes
(Newtonian) fluid (see Fung [29, Chapter 3]); however, in a vessel whose
characteristic dimension (diameter, e.g.) is about the same size as the
characteristic size of blood cells, blood behaves as a non-Newtonian fluid,
exhibiting shear-thinning and stress relaxation. Thurston [30, 31] pointed
out the viscoelastic behavior of blood while stating that the stress relaxation
is more significant for cases where the shear rate is low. With respect to
modeling blood, it has been observed that as the geometry of the flow changes,
the rheological characteristics of blood also change and as pointed out by
Humphrey and Delange [32, page 357], when blood flows in capillaries (5–8 μm in diameter), “the red blood cells go through
one at a time, with plasma in between.” Flow in larger vessels, in the range of
500–50 μm (see Fung [29, page 173]), exhibits a unique
phenomenon known as the Fahraeus-Lindqvist effect, whereby red blood cells
migrate towards the centerline, creating a thin layer of cell-depleted plasma
near the vessel wall. As they indicate, in such a case, the blood should be
treated as a solid and fluid mixture.In this paper,
we first provide a brief review of mixture theory, and then discuss certain
issues in constitutive modeling of blood. In the present formulation, we assume
blood to form a mixture consisting of RBCs suspended in plasma, while ignoring the
platelets, the white blood cells (WBCs), and the proteins in the sample. No
biochemical effects or interconversion of mass are considered in this model. The
volume fraction (or the concentration of the RBCs) is treated as a field
variable in this model. We further assume that the plasma behaves as a linearly
viscous fluid and the RBCs as an anisotropic nonlinear density-gradient-type
fluid. We obtain a constitutive relation for blood, based on the simplified
constitutive relations derived for plasma and RBCs. It is noted that we have only
discussed the development a model and that specific boundary value problems
need to be solved to test the efficacy of the model. A simple shear flow is studied
and an exact solution is obtained for a very special case; for more general
cases, it is necessary to solve the nonlinear coupled equations numerically.


## 2. A Brief Review of (a Two-component) Mixture Theory

In this section, we will  provide  a brief review of the
essential ideas and equations of a two-component mixture. At each instant of
time,t, it is assumed that each point in space is occupied by particles
belonging to both S1 and S2. Let X1 and X2 denote the
positions of particles of S1 and S2 in the reference
configuration, where (see Bowen [19], Truesdell [15])(2.1)x1=χ˜1(X1,t),x2=χ˜2(X2,t). These motions are assumed to be one-to-one, continuous, and
invertible. The kinematical quantities associated with these motions are(2.2)v1=D1χ˜1Dt,v2=D2χ˜2Dt,a1=D1v1Dt,a2=D2v2Dt,L1=∂v1∂x1,L2=∂v2∂x2,W1=12(L1−L1T),W2=12(L2−L2T),D1=12(L1+L1T),D2=12(L2+L2T), where v denotes
the velocity vector, a  the acceleration vector, L  the velocity gradient, D denotes the symmetric part of the
velocity gradient, and W the spin
tensor. D1/Dt denotes differentiation with respect to t, holding X1 fixed, and D2/Dt
denotes the same operation holding X2 fixed; in general for any scalar β, D(α)β/Dt=∂β/∂t+v(α)⋅gradβ, α=1,2,
and (for any vector w), D(α)w/Dt=∂w/∂t+(gradw)v(α).
Also, ρ1 and ρ2 are the bulk densities of the mixture components given by(2.3)ρ1=γρ10,ρ2=ϕρ20, where ρ10=ρf is the density of the first component (e.g.,
a fluid) in the reference configuration, ρ20=ρs is the density of the second component (e.g., solid
particles) in its reference configuration, γ is the volume fraction of the fluid component,
and ϕ is the volume fraction of the solid. For a saturated mixture γ=1−ϕ. The mixture density ρm is given by(2.4)ρm=ρ1+ρ2, and the mean velocity wm of the mixture is defined by(2.5)ρmwm=ρ1v1+ρ2v2. Once the individual stress tensors are derived (or proposed),
a mixture stress tensor can be defined as(2.6)Tm=T1+T2, where(2.7)T1=(1−ϕ)Tf,T2=Ts, so that the mixture stress tensor reduces to that of a pure
fluid as ϕ→0 and
to that of the solid particles (e.g., a densely packed granular material) as γ→0. T2 may also be written as T2=ϕT^s,
where T^s may be thought of as representing the stress
tensor in the reference configuration of the granular material. In the absence
of any thermochemical and electro-magnetic effects, the governing equations are
those of the conservation of mass and linear momentum.


### 2.1. Conservation of Mass

Assuming no interconversion of mass between the two
constituents, the equations for the conservation of mass for the two components
are(2.8)∂ρ1∂t+div(ρ1v1)=0,(2.9)∂ρ2∂t+div(ρ2v2)=0.


### 2.2. Conservation of Linear Momentum

LetT1 and T2 denote the partial
stress tensors. Then, the balance of
linear momentum equations for the two components is given by(2.10)ρ1D1v1Dt=divT1+ρ1b1+fI,(2.11)ρ2D2v2Dt=divT2+ρ2b2−fI, where b represents
the body force and fI represents the mechanical interaction (local exchange of momentum) between the
components.


### 2.3. Conservation of Angular
Momentum

The balance of moment of momentum implies that(2.12)T1+T2=T1T+T2T. The partial stresses however do not need to be symmetric. Equations
(2.8)–(2.11) represent
the basic governing equations for two-component flows, where the effects of electromagnetic
and temperature fields are ignored. These equations have to be supplemented
with constitutive relations for T1, T2, and fI. This is known as the
“closure” problem. In addition to the above balance equations (note that we
have ignored the balance of energy and the entropy equation), we need to
discuss certain constraints which might have to be considered in mixture theory.


### 2.4. Volume Additivity Constraint

Mills [33] derived the volume additivity
constraint by defining an incompressible mixture as one:(2.13)1ρ10dm1dV1+1ρ20dm2dV2=1, where “m”
denotes mass and “V” the volume, and where dV1+dV2=dV, dm1+dm2=dm, ρ10=dm1/dV1,
and ρ20=dm2/dV2.
Now, the densities in the current configuration, that is, in the mixture are
given by ρ1=dm1/dV, ρ2=dm2/dV.
Thus,(2.14)ρ1ρ10+ρ2ρ20=1. Other types of constraints such as
incompressibility and inextensibility can also be incorporated in mixture theory.


### 2.5. Surface (Boundary) Conditions

In continuum mechanics, in addition to the balance
(governing) equations, possibly subjected to some constraints, one must also
specify meaningful (physical) boundary conditions, in order to have a
well-posed problem. If the equations are nonlinear, the multiple solutions and
stability of those solutions are very often of interest. The typical boundary
conditions used in analytical/numerical procedures to solve any differential
equations are (i) Dirichlet boundary conditions
where the value of the unknowns is prescribed  on the boundary, (ii) Neumann
condition where the normal gradient of the unknowns is specified, and (iii)
boundary conditions where a combination of the unknown quantities and their
normal gradients is specified. The need for additional boundary conditions
arises in many areas of mechanics, whenever nonlinear or microstructural
theories are used. For example, in higher grade fluids in the mechanics
of non-Newtonian fluids (see Rajagopal and Kaloni [34]), in polar fluids (see
Atkin et al. [35]), in liquid crystals (see Leslie [36] and Ericksen
[37]), or in mixture theory (see Rajagopal et al. [38]), the need is
discussed and suggestions are offered.One of the fundamental questions in mixture theory is concerned
with the boundary conditions and how to split the (total) traction vector,
related to the (total) stress tensor, or the (total) velocity vector (see Ramtani [39]). In an
important paper, Rajagopal et al. [38] developed a novel scheme to split the
total stress for a class of problems in which the boundary of the mixture is
assumed to be in a state of saturation. These problems deal with diffusion of
fluids through nonlinear elastic materials, such as rubber. This additional
boundary condition is obtained by prescribing the variation of the Gibbs free
energy of dilution being zero. As a result of this thermodynamic restriction, a
relationship between the total stress tensor, the stretch tensor, and the
volume fraction of the solid component is obtained. This saturation condition
is explained by Rajagopal and Tao [3, page 31] as “…a state in which a small
element of the solid adjacent to the fluid is in a state in which it cannot
absorb any more fluid, that is whatever fluid enters the elemental volume along
the boundary has to exit through the elemental volume so that there is no
accumulation of the fluid.” Later, Tao and Rajagopal [3] suggested a purely
mechanical method for the splitting of the traction. If x is a point on the boundary surface of the mixture region and ts and tf are the surface tractions at x associated with the solid and the fluid components respectively,
then(2.15)ts=(Ts)Tn,tf=(Tf)Tn, where n is the
unit outward normal vector at x.
They assumed that the boundary surface fraction δSf/δS occupied by the fluid component at x is given by δSf/δS=δVf/δV, where δVf is the volume occupied by the fluid component
in δV.
Then, the traction vector tf at x is given by(2.16a)tf=δSfδSt, where t is the
total applied surface traction at x.
They also showed that(2.16b)(Tf)Tn=ρfρRft,(2.16c)(Ts)Tn=ρsρRst, where ρRs and ρRf are the mass densities of the solid component
and the fluid component in their reference states, and ρs and ρf denote their mass densities in the current
state. In Rajagopal and Massoudi’s model (see Massoudi et al. [12] for
details), it is not the tractions which are of interest in general but the
velocities. In their approach, the solid particles and fluid component are
assumed to both move and diffuse through each other, thereby transforming the
issue of traction boundary conditions to one of velocities at the boundary. For
free surface flows, on the other hand, the splitting of the traction vector
remains the main difficulty (see Ravindran et al. [40]). Sometimes these
additional boundary conditions can be provided from experimental data, and
sometimes they can be based on other theories or physical insights.In Rajagopal and Massoudi’s approach no couple stresses are
allowed. Nevertheless, due to the higher order gradients of volume fraction,
they also find it necessary to provide additional boundary conditions for
solving practical and simple boundary value problems (see Massoudi [41] for a
discussion of boundary conditions). For most practical applications, these can
sometimes be satisfied by certain symmetry conditions; in certain cases, the
values of the unknowns or their derivatives have to be specified as surface
conditions at the solid walls or at the free surface.For problems with bounded domains and with a certain degree
of symmetry, such as fully developed flow in a vertical pipe (see Gudhe et al.
[42]), one can prescribe at the center of the pipe the following condition:(2.17a)dϕdr=0,(2.17b)dvdr=dudr=0, where u and v are the velocities. It is much more difficult
to specify the value of the volume fraction (or concentration) at the solid
surfaces. For example, based on experimental evidence, or artificially and
manually gluing particles (or a thin layer) to the wall, one can prescribe a
value such as(2.18a)ϕ=ϕw(at the wall). Other helpful conditions, though strictly speaking not a
boundary condition, are integral conditions such as a mass flux, which provide
an average value. For example, for the case of pipe flow, we can use(2.18b)N=2π∫0rϕrdr, where the value of N is an input to the problem. In most
cases, it can be assumed that the no-slip boundary condition applies to the
velocity components. Therefore,(2.19a)u=v=0at the boundaries. However, in many situations slip may
occur at the wall, and therefore the classical assumption of adherence boundary
condition at the wall no longer applies. For most applications, a
generalization of Navier’s hypothesis can be used (see Rajagopal [43])(2.19b)u⋅t=−k(Tn⋅t),k>0, where u is the velocity vector of the fluid, T the stress tensor for the fluid, t and n are the unit tangent and normal vectors at the boundary, and k,
in general, can depend on the normal stress as well as the shear rate. Massoudi
and Phuoc [44] proposed that for granular materials the slip velocity is
proportional to the stress vector at the wall, that is, us=g[(Tsn)x,(Tsn)y], where Ts is the stress
tensor for the granular component, n is the unit normal vector, and g in general could be a function of surface roughness, volume fraction (density),
shear rate, and so forth.For free surface flows, in general, the location of the free
surface is not known and one must find this surface as part of the solution.
For steady flows, one can use the kinematical constraint(2.20)u⋅n=0. For unsteady flows, the problem of specifying the boundary
conditions is more complicated. On the free surface it is known that the stress
is zero, that is, one has the traction-free boundary condition. For a mixture,
however, it is not clear whether the total traction vector t is zero, or whether the individual traction vectors ts and tf are zero. For a fully developed flow, where the
location of the free surface, that is, the constant height of the free surface
is known a priori, and based on
the work of Tao and Rajagopal [45], Ravindran et al. [40] assumed that the
tangential components of the individual traction vectors are zero, while the
normal components are weighted according to the volume fraction. That is, they
assumed that the atmospheric pressure on the top surface is split between the
two components in the ratio of their respective volume fractions:(2.21a)Ts(01)=(0ϕPatm), and for the fluid, we have(2.21b)Tf(01)=(0(1−ϕ)Patm). In general, if one attempts to eliminate the pressure term
via cross-differentiation, the order of the equations increases, and one needs additional
boundary condition. This can be obtained by specifying the flow rate of the
mixture (see Beevers and Craine [46]):(2.22a)Qvolumetric=∫01VmdY=∫01[(1−ϕ)u+ϕv]dY. Alternatively, the flow rate of the mixture can be
prescribed. The mass flow rate for a two-component mixture is defined as(2.22b)Qm=∫01ρmVmdY=∫01[(1−ϕ)ρf+ϕρs][(1−ϕ)u+ϕv]dY. When the mixture is neutrally buoyant, that is, ρf=ρs=ρ,
these two quantities are related by Qm=ρQ.Finally, for a complete study of a thermomechanical problem,
not only in mixture theory, but in continuum mechanics in general, the second law
of thermodynamics has to be considered. In other words, in addition to other
principles in continuum mechanics such as material symmetry, and frame
indifference, the second law imposes important  restrictions on the type of motion and/or the
constitutive parameters (for a thorough discussion of important concepts in constitutive
equations of mechanics, we refer the reader to the books by Antman [47]
Maugin [48], Coussot [49], Liu [50], Batra [51]). Since, there is no
general agreement on the functional form of the constitutive relation and since
the Helmholtz free energy is not known, a complete thermodynamical treatment of
the model used in our studies is lacking. In recent years, Rajagopal and
colleagues (see, e.g., Rajagopal and Srinivasa [52, 53]) have devised a
thermodynamic framework, the multiple natural configuration theory, by appealing to the maximization of the rate
of entropy production to obtain a class of constitutive relations for many
different types of materials. Unlike the traditional thermodynamic approach
whereby a form for the stress is assumed (or derived) and restriction on the
material parameters is obtained by invoking the Clausius-Duhem inequality, in
their thermodynamic framework, they assume specific forms for the Helmholtz
potential and the rate of dissipation—reflecting on how
the energy is stored in the body and the way in which the body dissipates it. In
the next section, we discuss a general method on how to obtain the balance equations
of mixture if the balance equations of each component are given.


## 2.1. Conservation of Mass

Assuming no interconversion of mass between the two
constituents, the equations for the conservation of mass for the two components
are(2.8)∂ρ1∂t+div(ρ1v1)=0,(2.9)∂ρ2∂t+div(ρ2v2)=0.


## 2.2. Conservation of Linear Momentum

LetT1 and T2 denote the partial
stress tensors. Then, the balance of
linear momentum equations for the two components is given by(2.10)ρ1D1v1Dt=divT1+ρ1b1+fI,(2.11)ρ2D2v2Dt=divT2+ρ2b2−fI, where b represents
the body force and fI represents the mechanical interaction (local exchange of momentum) between the
components.


## 2.3. Conservation of Angular
Momentum

The balance of moment of momentum implies that(2.12)T1+T2=T1T+T2T. The partial stresses however do not need to be symmetric. Equations
(2.8)–(2.11) represent
the basic governing equations for two-component flows, where the effects of electromagnetic
and temperature fields are ignored. These equations have to be supplemented
with constitutive relations for T1, T2, and fI. This is known as the
“closure” problem. In addition to the above balance equations (note that we
have ignored the balance of energy and the entropy equation), we need to
discuss certain constraints which might have to be considered in mixture theory.


## 2.4. Volume Additivity Constraint

Mills [33] derived the volume additivity
constraint by defining an incompressible mixture as one:(2.13)1ρ10dm1dV1+1ρ20dm2dV2=1, where “m”
denotes mass and “V” the volume, and where dV1+dV2=dV, dm1+dm2=dm, ρ10=dm1/dV1,
and ρ20=dm2/dV2.
Now, the densities in the current configuration, that is, in the mixture are
given by ρ1=dm1/dV, ρ2=dm2/dV.
Thus,(2.14)ρ1ρ10+ρ2ρ20=1. Other types of constraints such as
incompressibility and inextensibility can also be incorporated in mixture theory.


## 2.5. Surface (Boundary) Conditions

In continuum mechanics, in addition to the balance
(governing) equations, possibly subjected to some constraints, one must also
specify meaningful (physical) boundary conditions, in order to have a
well-posed problem. If the equations are nonlinear, the multiple solutions and
stability of those solutions are very often of interest. The typical boundary
conditions used in analytical/numerical procedures to solve any differential
equations are (i) Dirichlet boundary conditions
where the value of the unknowns is prescribed  on the boundary, (ii) Neumann
condition where the normal gradient of the unknowns is specified, and (iii)
boundary conditions where a combination of the unknown quantities and their
normal gradients is specified. The need for additional boundary conditions
arises in many areas of mechanics, whenever nonlinear or microstructural
theories are used. For example, in higher grade fluids in the mechanics
of non-Newtonian fluids (see Rajagopal and Kaloni [34]), in polar fluids (see
Atkin et al. [35]), in liquid crystals (see Leslie [36] and Ericksen
[37]), or in mixture theory (see Rajagopal et al. [38]), the need is
discussed and suggestions are offered.One of the fundamental questions in mixture theory is concerned
with the boundary conditions and how to split the (total) traction vector,
related to the (total) stress tensor, or the (total) velocity vector (see Ramtani [39]). In an
important paper, Rajagopal et al. [38] developed a novel scheme to split the
total stress for a class of problems in which the boundary of the mixture is
assumed to be in a state of saturation. These problems deal with diffusion of
fluids through nonlinear elastic materials, such as rubber. This additional
boundary condition is obtained by prescribing the variation of the Gibbs free
energy of dilution being zero. As a result of this thermodynamic restriction, a
relationship between the total stress tensor, the stretch tensor, and the
volume fraction of the solid component is obtained. This saturation condition
is explained by Rajagopal and Tao [3, page 31] as “…a state in which a small
element of the solid adjacent to the fluid is in a state in which it cannot
absorb any more fluid, that is whatever fluid enters the elemental volume along
the boundary has to exit through the elemental volume so that there is no
accumulation of the fluid.” Later, Tao and Rajagopal [3] suggested a purely
mechanical method for the splitting of the traction. If x is a point on the boundary surface of the mixture region and ts and tf are the surface tractions at x associated with the solid and the fluid components respectively,
then(2.15)ts=(Ts)Tn,tf=(Tf)Tn, where n is the
unit outward normal vector at x.
They assumed that the boundary surface fraction δSf/δS occupied by the fluid component at x is given by δSf/δS=δVf/δV, where δVf is the volume occupied by the fluid component
in δV.
Then, the traction vector tf at x is given by(2.16a)tf=δSfδSt, where t is the
total applied surface traction at x.
They also showed that(2.16b)(Tf)Tn=ρfρRft,(2.16c)(Ts)Tn=ρsρRst, where ρRs and ρRf are the mass densities of the solid component
and the fluid component in their reference states, and ρs and ρf denote their mass densities in the current
state. In Rajagopal and Massoudi’s model (see Massoudi et al. [12] for
details), it is not the tractions which are of interest in general but the
velocities. In their approach, the solid particles and fluid component are
assumed to both move and diffuse through each other, thereby transforming the
issue of traction boundary conditions to one of velocities at the boundary. For
free surface flows, on the other hand, the splitting of the traction vector
remains the main difficulty (see Ravindran et al. [40]). Sometimes these
additional boundary conditions can be provided from experimental data, and
sometimes they can be based on other theories or physical insights.In Rajagopal and Massoudi’s approach no couple stresses are
allowed. Nevertheless, due to the higher order gradients of volume fraction,
they also find it necessary to provide additional boundary conditions for
solving practical and simple boundary value problems (see Massoudi [41] for a
discussion of boundary conditions). For most practical applications, these can
sometimes be satisfied by certain symmetry conditions; in certain cases, the
values of the unknowns or their derivatives have to be specified as surface
conditions at the solid walls or at the free surface.For problems with bounded domains and with a certain degree
of symmetry, such as fully developed flow in a vertical pipe (see Gudhe et al.
[42]), one can prescribe at the center of the pipe the following condition:(2.17a)dϕdr=0,(2.17b)dvdr=dudr=0, where u and v are the velocities. It is much more difficult
to specify the value of the volume fraction (or concentration) at the solid
surfaces. For example, based on experimental evidence, or artificially and
manually gluing particles (or a thin layer) to the wall, one can prescribe a
value such as(2.18a)ϕ=ϕw(at the wall). Other helpful conditions, though strictly speaking not a
boundary condition, are integral conditions such as a mass flux, which provide
an average value. For example, for the case of pipe flow, we can use(2.18b)N=2π∫0rϕrdr, where the value of N is an input to the problem. In most
cases, it can be assumed that the no-slip boundary condition applies to the
velocity components. Therefore,(2.19a)u=v=0at the boundaries. However, in many situations slip may
occur at the wall, and therefore the classical assumption of adherence boundary
condition at the wall no longer applies. For most applications, a
generalization of Navier’s hypothesis can be used (see Rajagopal [43])(2.19b)u⋅t=−k(Tn⋅t),k>0, where u is the velocity vector of the fluid, T the stress tensor for the fluid, t and n are the unit tangent and normal vectors at the boundary, and k,
in general, can depend on the normal stress as well as the shear rate. Massoudi
and Phuoc [44] proposed that for granular materials the slip velocity is
proportional to the stress vector at the wall, that is, us=g[(Tsn)x,(Tsn)y], where Ts is the stress
tensor for the granular component, n is the unit normal vector, and g in general could be a function of surface roughness, volume fraction (density),
shear rate, and so forth.For free surface flows, in general, the location of the free
surface is not known and one must find this surface as part of the solution.
For steady flows, one can use the kinematical constraint(2.20)u⋅n=0. For unsteady flows, the problem of specifying the boundary
conditions is more complicated. On the free surface it is known that the stress
is zero, that is, one has the traction-free boundary condition. For a mixture,
however, it is not clear whether the total traction vector t is zero, or whether the individual traction vectors ts and tf are zero. For a fully developed flow, where the
location of the free surface, that is, the constant height of the free surface
is known a priori, and based on
the work of Tao and Rajagopal [45], Ravindran et al. [40] assumed that the
tangential components of the individual traction vectors are zero, while the
normal components are weighted according to the volume fraction. That is, they
assumed that the atmospheric pressure on the top surface is split between the
two components in the ratio of their respective volume fractions:(2.21a)Ts(01)=(0ϕPatm), and for the fluid, we have(2.21b)Tf(01)=(0(1−ϕ)Patm). In general, if one attempts to eliminate the pressure term
via cross-differentiation, the order of the equations increases, and one needs additional
boundary condition. This can be obtained by specifying the flow rate of the
mixture (see Beevers and Craine [46]):(2.22a)Qvolumetric=∫01VmdY=∫01[(1−ϕ)u+ϕv]dY. Alternatively, the flow rate of the mixture can be
prescribed. The mass flow rate for a two-component mixture is defined as(2.22b)Qm=∫01ρmVmdY=∫01[(1−ϕ)ρf+ϕρs][(1−ϕ)u+ϕv]dY. When the mixture is neutrally buoyant, that is, ρf=ρs=ρ,
these two quantities are related by Qm=ρQ.Finally, for a complete study of a thermomechanical problem,
not only in mixture theory, but in continuum mechanics in general, the second law
of thermodynamics has to be considered. In other words, in addition to other
principles in continuum mechanics such as material symmetry, and frame
indifference, the second law imposes important  restrictions on the type of motion and/or the
constitutive parameters (for a thorough discussion of important concepts in constitutive
equations of mechanics, we refer the reader to the books by Antman [47]
Maugin [48], Coussot [49], Liu [50], Batra [51]). Since, there is no
general agreement on the functional form of the constitutive relation and since
the Helmholtz free energy is not known, a complete thermodynamical treatment of
the model used in our studies is lacking. In recent years, Rajagopal and
colleagues (see, e.g., Rajagopal and Srinivasa [52, 53]) have devised a
thermodynamic framework, the multiple natural configuration theory, by appealing to the maximization of the rate
of entropy production to obtain a class of constitutive relations for many
different types of materials. Unlike the traditional thermodynamic approach
whereby a form for the stress is assumed (or derived) and restriction on the
material parameters is obtained by invoking the Clausius-Duhem inequality, in
their thermodynamic framework, they assume specific forms for the Helmholtz
potential and the rate of dissipation—reflecting on how
the energy is stored in the body and the way in which the body dissipates it. In
the next section, we discuss a general method on how to obtain the balance equations
of mixture if the balance equations of each component are given.


## 3. Equations of Motion for The Mixture

To obtain the
conservation of mass for the mixture, we add (2.8) and (2.9) to obtain(3.1)∂(ρ1+ρ2)∂t+div(ρ1v1+ρ2v2)=0. Now, the
mass-weighted velocity of the mixture is defined as (see Bowen [19])(3.2)ρmwm=ρ1v1+ρ2v2. Thus, substituting
this equation and recalling that ρm=ρ1+ρ2,
we obtain the conservation of mass for the mixture:(3.3)∂ρm∂t+div(ρmwm)=0. Similarly, to obtain
the conservation of linear momentum for the mixture, we can add (2.10) and (2.11):(3.4)ρ1D1v1Dt+ρ2D2v2Dt=div(T1+T2)+ρ1b1+ρ2b2. Defining(3.5)ρmbm=ρ1b1+ρ2b2. Equation (3.4) is
rewritten as, using (2.6),(3.6)ρ1D1v1Dt+ρ2D2v2Dt=divTm+ρmbm. This is the balance of linear momentum for
the mixture, and unlike (3.3) which looks similar to the conservation of mass
for a single continuum, this equation appears different from that of the linear
momentum equation for a single component continuum, due to the two time
derivatives on the left-hand side of the equation. An alternative form for the
balance of linear momentum for the mixture as a whole, as prescribed in the
classical mixture theory (see, e.g., Müller [54], Bowen [19]), is given by(3.7)ρmDwmDt=divTmix+ρmbm, where now(3.8)Tmix=T(1)+T(2)−∑α=12ραq(α)⊗q(α), where q(α) is the diffusion velocity vector for the αth constituent at (x,t) defined by(3.9)q(α)=V(α)−w, where(3.10a)V(1)=u,(3.10b)V(2)=v. Thus,(3.11)Tmix=T(1)+T(2)−[ρ1(u−w)⊗(u−w)+ρ2(v−w)⊗(v−w)]. It is clearly seen that with this definition
for mixture stress tensor, the balance of linear momentum for the mixture,
(3.7), would appear differently, since the terms included in the bracket on the
right-hand side of (3.11) now would appear in (3.7). As indicated by Massoudi [16],
this form is not very amenable to practical engineering problems where one is
interested in mixture properties, such as viscosity. As mentioned earlier, another
way of defining the mixture stress tensor is due to Green and Naghdi [55, 56], where(3.12)Tm=T(1)+T(2),where T(1) and T(2) are the constituent (partial) stress tensor
relations for the two components as used  in the above equations. Whether we use
(3.11) or (3.12), the specific forms of T(1) and T(2) depend upon the way in which the constitutive
relations are obtained and the types of restrictions imposed on the materials
properties. The advantage of presenting the balance equations for the mixture,
(3.3) and (3.4), as opposed to the balance equations for each component, is that
in this formulation the interaction forces cancel each other and one no longer
needs to consider constitutive relations for interaction forces. The disadvantage
of this approach is that one can no longer obtain detailed information about
velocity, pressure, and concentration distributions.In the next section,
we will discuss the approach of Massoudi and Rajagopal, as one of the possible
approaches within mixture theory, in formulating constitutive equations for the
stress tensors and the interaction forces.


## 4. Constitutive Relations

Mathematically, the purpose of the constitutive relations is
to supply connections between kinematical, mechanical, electromagnetic, and
thermal fields that are compatible with the balance equations and that, in
conjunction with them, provide a  theory that is solvable for properly posed
problems. Deriving constitutive relations for the stress tensors and the
interaction forces is among
the outstanding issues of research in multicomponent flows. In general, the
constitutive expressions forT1 and T2 depend on the
kinematical quantities associated with both the constituents. However, it can be assumed that T1 and T2 depend only on the kinematical quantities associated
with the particles (component 1)  and fluid
(component 2), respectively (sometimes called the principle of component
separation, see Adkins [57, 58]).In this part of the paper, we provide a brief description of
Massoudi and Rajagopal [M-R] approach. We then modify the M-R model to obtain a
constitutive relation for blood. The M-R model considers a mixture of an
incompressible fluid infused with solid particles, wherein the principles of mechanics
of granular materials are used to describe the behavior of particles. This
model has been used and discussed in Massoudi [2, 16, 41, 59–62], Johnson et al. [9, 10], Rajagopal et al. [63–65], Massoudi
et al. [12], Massoudi and Johnson [66], Massoudi and Lakshmana Rao [67], and
Ravindran et al. [40].


### 4.1. Stress Tensor for the Fluid Component

In most practical engineering problems studied by Massoudi
and Rajagopal, the fluid can be assumed to behave as a linearly viscous fluid:(4.1)T1=[−p(ρ1)+λ1(ρ1)trD1]1+2μ1(ρ1)D1, where p is the fluid pressure, μ1 is the viscosity, and D1 is the symmetric part of the velocity gradient
of the fluid, and λ1 is the second coefficient of viscosity. Obviously,
this assumption can be modified if the fluid component is a nonlinear
(non-Newtonian) fluid.


### 4.2. Stress Tensor for the Solid Component

The stress tensorT2  in a flowing granular material may depend on
the manner in which the granular material is distributed, that is, the volume fraction ϕ and possibly also its gradient, and the
symmetric part of the velocity gradient tensor D2. Based on this observation,
they assume that (see also Cowin [68])(4.2)T2=f(ϕ,∇ϕ,D2). Using
standard arguments in mechanics, restrictions can be found on the form of the
above constitutive expression based on the assumption of frame-indifference,
isotropy, and so forth. There
could be further restrictions on the form of the constitutive expression because
of internal constraints, such as, incompressibility and thermodynamics
restrictions due to Clausius-Duhem inequality. A constitutive model that
predicts the possibility of normal stress-differences and is properly frame
invariant was given by Rajagopal and Massoudi [69](4.3)T2=[βo+β1∇ϕ⋅∇ϕ+β2trD2]1+β3D2+β4∇ϕ⊗∇ϕ+β5D22, where(4.4)D2=12[∇v2+(∇v2)T], where ρ2=ρsϕ,
with ρs being constant, and the β's are material properties, which in general are
functions of the appropriate principal invariants of the density gradient and
the symmetric part of the velocity gradient. They assumed the following forms:(4.5)βo=kϕ,k<0,(4.6)β1=β1*(1+ϕ+ϕ2),β2=β2*(ϕ+ϕ2),β3=β3*(ϕ+ϕ2),β4=β4*(1+ϕ+ϕ2),β5=β5*(1+ϕ+ϕ2). The above representation can be viewed as Taylor
series approximation for the material
parameters (Rajagopal et al. [65]). Such a quadratic dependence, at least for the viscosity β3,
is on the basis of dynamic simulations of particle interactions (Walton and
Braun [70, 71]). Further restrictions
on the coefficients have been obtained by using the argument that the stress
should vanish as ϕ→0, therefore(4.7)β50=β30=β20. This is a principle of the limiting case, that is, if there
are no particles, then ϕ and gradϕ are zero, and the stress should be zero. Since,
the kinematical terms D2, D22
, and tr D2 multiplied by β2, β3, and β5 do
not necessarily go to zero when there are no particles, the above restrictions
are necessary. Furthermore, Rajagopal
and Massoudi [69] and Rajagopal et al. [63] have shown that(4.8)β1+β4>0,k<0, as compression should lead to densification of the material.
They suggested the following rheological interpretation to the material
parameters: β0(ϕ) is similar to pressure in a compressible fluid
and is to be given by an equation of state, β2(ϕ) is like the second coefficient of viscosity in
a compressible fluid, β1(ϕ) and β4(ϕ) are the material parameters connected with the
distribution of the granular materials, β3(ϕ) is the viscosity of the granular materials,
and β5(ϕ) is similar to the viscosity term in the
Reiner-Rivlin model (Reiner [72, 73], Rivlin [74]), often referred to as the cross-viscosity. For the rheological
parameters, β2, β3,
and β5,
one can use the methods available in the
mechanics of non-Newtonian fluids to find out more information about the signs.
Obviously, since β3 is related to the shear viscosity, it is
assumed to be positive. The simulation studies of Walton and Braun [70, 71]
suggest a structure of β3 similar to the one assumed by Rajagopal and
Massoudi. As shown in Rajagopal et al. [65] by using an orthogonal rheometer,
and measuring the forces and moments exerted on the disks, one can characterize
the materials moduli β's.
Rajagopal et al. [75] and Baek et al. [76] discuss the details of
experimental techniques using orthogonal and torsional rheometers to measure
the material properties β1 and β4.
In a certain sense, the model represented by (4.3) has the same structure as the
Reiner-Rivlin model, whereby the effects of density (or volume fraction)
gradients are also included. This model has been used in a variety of simple
applications such as inclined flow and flow in a vertical pipe.Looking at (4.3) with (4.5)–(4.7), it can be
shown that this model is capable of predicting both normal stress differences
in a simple shear flow problem (for details see Massoudi and Mehrabadi [77]).
The volume fraction field ϕ(x,t) plays a major role in this model. That is, even though we talk of distinct
solid particles with a certain diameter, in this theory, the particles through
the introduction of the volume fraction field have been homogenized. In other
words, it is assumed that the material properties of the ensemble are
continuous functions of position. That
is, the material may be divided indefinitely without losing any of its defining
properties (see also Collins [78]). A distributed volume Vt=∫ϕdV and a distributed mass M=∫ρsϕdV can be defined, where the function ϕ is an independent kinematical variable called
the volume distribution function and has the property(4.9)0≤ϕ(x,t)<ϕmax<1. The function ϕ is represented as a continuous function of
position and time; in reality, ϕ in a granular system is either one or zero at
any position and time, depending upon whether one is pointing to a granule or
to the void space at that position. That
is, the real volume distribution content has been averaged, in some sense, over
the neighborhood of any given position. It should be mentioned that in practice ϕ is never equal to one; its maximum value,
generally designated as the maximum packing fraction, depends on the shape,
size, method of packing, and so forth. (For a review of constitutive modelling
of (flowing) granular materials, see Massoudi [79]). This further complicates
experimental identification of constitutive parameters, since it is not
possible to study a “pure” sample of the granular phase.


### 4.3. Constitutive Relation for the Interaction Forces

The interaction force, in general, depends on the fluid
pressure gradient, the density gradient (the buoyancy forces), the relative
velocity (the drag force on the particles), the relative acceleration (the
virtual mass of the particles), the magnitude of the rate of the deformation
tensor of the fluid (the lift force on the particles), the spinning motion, as
well as the translation of particles (the Faxen’s force), the particles
tendency to move toward the region of higher velocity (the Magnus force), the
history of the particle motion (the Bassett force), the temperature gradient,
and so forth. For laminar flow of a mixture of an incompressible fluid infused
with particles, the mechanical interaction force can be assumed to be of the
form (Johnson et al. [80])(4.10)fI=A1gradϕ+A2F(ϕ)(v2-v1)+A3ϕ(2trD12)−1/4D1(v2-v1)+A4ϕ(W2-W1)(v2-v1)+A5avm, where W2 and W1 are the spin tensors for the two components.
The terms on the right-hand side of this equation reflect the presence of (nonuniform)
concentration distribution (diffusion), drag, slip-shear lift, spin lift, and
virtual mass. If the flow is assumed to be steady, the virtual mass effects can
be neglected. Most of these coefficients have not been measured or extensively
studied for general two-component flows. Johnson et al. [80], based on
previous experimental observations and certain analytical approximations,
suggested the following coefficients:(4.11)A2=92μfa2;A3=3(6.46)4πρf1/2μf1/2a,A4=34ρf,A5=2π3a3ϕ1+2ϕ1-ϕ. Most of these coefficients are given either as
a generalization of single particle result or as a result of some other
limiting conditions. If the particles are nonspherical, such as fibers, or a
blood cell then the directionality may become an important element, and in that
case we need to use (or develop) the constitutive relations which take the
microstructure and directionality, that is, anisotropic nature of the material,
into account (cf., Ericksen [37, 81], Leslie [82]). There are obvious
“gray” areas that await further studies for clarification. For example, there
have been no investigations as to what form A1 may have. The remaining coefficients have not
been extensively studied for general two-component flows; thus, the forms given
above are ad hoc applications
of results that are strictly valid under more restricted conditions. Despite
the assumptions involved, these expressions effectively describe how the
interaction forces vary with the system parameters. The coefficients can be
considered functions of ϕ only (because ϕ varies with position in the flow) for the
purpose of performing numerical studies. The way in which Massoudi and
Rajagopal use (4.10) along with (4.11) is that since the governing equations are
made dimensionless, the material parameters such as viscosity and density along
with the numerical coefficients appearing in these correlations are absorbed in
the dimensionless numbers. Therefore, the actual numerical values of these
experimentally obtained coefficients do not enter into their calculations
directly; instead they perform a parametric study for a range of these
dimensionless numbers. Massoudi [2, 62] has recently reviewed the subject
of interaction mechanisms in multicomponent flows. In the next section, we
discuss a possible extension of the fluid-solid mixture theory formulation of
Massoudi and Rajagopal to blood.


## 4.1. Stress Tensor for the Fluid Component

In most practical engineering problems studied by Massoudi
and Rajagopal, the fluid can be assumed to behave as a linearly viscous fluid:(4.1)T1=[−p(ρ1)+λ1(ρ1)trD1]1+2μ1(ρ1)D1, where p is the fluid pressure, μ1 is the viscosity, and D1 is the symmetric part of the velocity gradient
of the fluid, and λ1 is the second coefficient of viscosity. Obviously,
this assumption can be modified if the fluid component is a nonlinear
(non-Newtonian) fluid.


## 4.2. Stress Tensor for the Solid Component

The stress tensorT2  in a flowing granular material may depend on
the manner in which the granular material is distributed, that is, the volume fraction ϕ and possibly also its gradient, and the
symmetric part of the velocity gradient tensor D2. Based on this observation,
they assume that (see also Cowin [68])(4.2)T2=f(ϕ,∇ϕ,D2). Using
standard arguments in mechanics, restrictions can be found on the form of the
above constitutive expression based on the assumption of frame-indifference,
isotropy, and so forth. There
could be further restrictions on the form of the constitutive expression because
of internal constraints, such as, incompressibility and thermodynamics
restrictions due to Clausius-Duhem inequality. A constitutive model that
predicts the possibility of normal stress-differences and is properly frame
invariant was given by Rajagopal and Massoudi [69](4.3)T2=[βo+β1∇ϕ⋅∇ϕ+β2trD2]1+β3D2+β4∇ϕ⊗∇ϕ+β5D22, where(4.4)D2=12[∇v2+(∇v2)T], where ρ2=ρsϕ,
with ρs being constant, and the β's are material properties, which in general are
functions of the appropriate principal invariants of the density gradient and
the symmetric part of the velocity gradient. They assumed the following forms:(4.5)βo=kϕ,k<0,(4.6)β1=β1*(1+ϕ+ϕ2),β2=β2*(ϕ+ϕ2),β3=β3*(ϕ+ϕ2),β4=β4*(1+ϕ+ϕ2),β5=β5*(1+ϕ+ϕ2). The above representation can be viewed as Taylor
series approximation for the material
parameters (Rajagopal et al. [65]). Such a quadratic dependence, at least for the viscosity β3,
is on the basis of dynamic simulations of particle interactions (Walton and
Braun [70, 71]). Further restrictions
on the coefficients have been obtained by using the argument that the stress
should vanish as ϕ→0, therefore(4.7)β50=β30=β20. This is a principle of the limiting case, that is, if there
are no particles, then ϕ and gradϕ are zero, and the stress should be zero. Since,
the kinematical terms D2, D22
, and tr D2 multiplied by β2, β3, and β5 do
not necessarily go to zero when there are no particles, the above restrictions
are necessary. Furthermore, Rajagopal
and Massoudi [69] and Rajagopal et al. [63] have shown that(4.8)β1+β4>0,k<0, as compression should lead to densification of the material.
They suggested the following rheological interpretation to the material
parameters: β0(ϕ) is similar to pressure in a compressible fluid
and is to be given by an equation of state, β2(ϕ) is like the second coefficient of viscosity in
a compressible fluid, β1(ϕ) and β4(ϕ) are the material parameters connected with the
distribution of the granular materials, β3(ϕ) is the viscosity of the granular materials,
and β5(ϕ) is similar to the viscosity term in the
Reiner-Rivlin model (Reiner [72, 73], Rivlin [74]), often referred to as the cross-viscosity. For the rheological
parameters, β2, β3,
and β5,
one can use the methods available in the
mechanics of non-Newtonian fluids to find out more information about the signs.
Obviously, since β3 is related to the shear viscosity, it is
assumed to be positive. The simulation studies of Walton and Braun [70, 71]
suggest a structure of β3 similar to the one assumed by Rajagopal and
Massoudi. As shown in Rajagopal et al. [65] by using an orthogonal rheometer,
and measuring the forces and moments exerted on the disks, one can characterize
the materials moduli β's.
Rajagopal et al. [75] and Baek et al. [76] discuss the details of
experimental techniques using orthogonal and torsional rheometers to measure
the material properties β1 and β4.
In a certain sense, the model represented by (4.3) has the same structure as the
Reiner-Rivlin model, whereby the effects of density (or volume fraction)
gradients are also included. This model has been used in a variety of simple
applications such as inclined flow and flow in a vertical pipe.Looking at (4.3) with (4.5)–(4.7), it can be
shown that this model is capable of predicting both normal stress differences
in a simple shear flow problem (for details see Massoudi and Mehrabadi [77]).
The volume fraction field ϕ(x,t) plays a major role in this model. That is, even though we talk of distinct
solid particles with a certain diameter, in this theory, the particles through
the introduction of the volume fraction field have been homogenized. In other
words, it is assumed that the material properties of the ensemble are
continuous functions of position. That
is, the material may be divided indefinitely without losing any of its defining
properties (see also Collins [78]). A distributed volume Vt=∫ϕdV and a distributed mass M=∫ρsϕdV can be defined, where the function ϕ is an independent kinematical variable called
the volume distribution function and has the property(4.9)0≤ϕ(x,t)<ϕmax<1. The function ϕ is represented as a continuous function of
position and time; in reality, ϕ in a granular system is either one or zero at
any position and time, depending upon whether one is pointing to a granule or
to the void space at that position. That
is, the real volume distribution content has been averaged, in some sense, over
the neighborhood of any given position. It should be mentioned that in practice ϕ is never equal to one; its maximum value,
generally designated as the maximum packing fraction, depends on the shape,
size, method of packing, and so forth. (For a review of constitutive modelling
of (flowing) granular materials, see Massoudi [79]). This further complicates
experimental identification of constitutive parameters, since it is not
possible to study a “pure” sample of the granular phase.


## 4.3. Constitutive Relation for the Interaction Forces

The interaction force, in general, depends on the fluid
pressure gradient, the density gradient (the buoyancy forces), the relative
velocity (the drag force on the particles), the relative acceleration (the
virtual mass of the particles), the magnitude of the rate of the deformation
tensor of the fluid (the lift force on the particles), the spinning motion, as
well as the translation of particles (the Faxen’s force), the particles
tendency to move toward the region of higher velocity (the Magnus force), the
history of the particle motion (the Bassett force), the temperature gradient,
and so forth. For laminar flow of a mixture of an incompressible fluid infused
with particles, the mechanical interaction force can be assumed to be of the
form (Johnson et al. [80])(4.10)fI=A1gradϕ+A2F(ϕ)(v2-v1)+A3ϕ(2trD12)−1/4D1(v2-v1)+A4ϕ(W2-W1)(v2-v1)+A5avm, where W2 and W1 are the spin tensors for the two components.
The terms on the right-hand side of this equation reflect the presence of (nonuniform)
concentration distribution (diffusion), drag, slip-shear lift, spin lift, and
virtual mass. If the flow is assumed to be steady, the virtual mass effects can
be neglected. Most of these coefficients have not been measured or extensively
studied for general two-component flows. Johnson et al. [80], based on
previous experimental observations and certain analytical approximations,
suggested the following coefficients:(4.11)A2=92μfa2;A3=3(6.46)4πρf1/2μf1/2a,A4=34ρf,A5=2π3a3ϕ1+2ϕ1-ϕ. Most of these coefficients are given either as
a generalization of single particle result or as a result of some other
limiting conditions. If the particles are nonspherical, such as fibers, or a
blood cell then the directionality may become an important element, and in that
case we need to use (or develop) the constitutive relations which take the
microstructure and directionality, that is, anisotropic nature of the material,
into account (cf., Ericksen [37, 81], Leslie [82]). There are obvious
“gray” areas that await further studies for clarification. For example, there
have been no investigations as to what form A1 may have. The remaining coefficients have not
been extensively studied for general two-component flows; thus, the forms given
above are ad hoc applications
of results that are strictly valid under more restricted conditions. Despite
the assumptions involved, these expressions effectively describe how the
interaction forces vary with the system parameters. The coefficients can be
considered functions of ϕ only (because ϕ varies with position in the flow) for the
purpose of performing numerical studies. The way in which Massoudi and
Rajagopal use (4.10) along with (4.11) is that since the governing equations are
made dimensionless, the material parameters such as viscosity and density along
with the numerical coefficients appearing in these correlations are absorbed in
the dimensionless numbers. Therefore, the actual numerical values of these
experimentally obtained coefficients do not enter into their calculations
directly; instead they perform a parametric study for a range of these
dimensionless numbers. Massoudi [2, 62] has recently reviewed the subject
of interaction mechanisms in multicomponent flows. In the next section, we
discuss a possible extension of the fluid-solid mixture theory formulation of
Massoudi and Rajagopal to blood.


## 5. Constitutive Relations for Plasma and Red Blood Cells

The
physical and biological processes governing the flow of blood are intimately responsible
for safety and efficacy of all blood-wetted medical devices. The quest to
design improved cardiovascular devices is
however stifled by the inadequacies of current understanding of blood trauma
and thrombosis. Contemporary design relies upon formula for blood describing
(1) rheology, (2) cell trauma (hemolysis), and (3) thrombosis, that are based
primarily on empiricism. These relations are application-specific at best, and
are moredescriptive than predictive. Furthermore, we now
understand that these three phenomena are more closely coupled than that previously
appreciated. Unfortunately, the deficiencies of accurate predictive
mathematical models have
prevented engineering rigor to supplant the common practices which rely on
intuition and experimental trial-and-error. Chief on the list of biological
phenomena that are difficult to predict a priori is thrombosis and hemolysis,
collectively referred to as blood trauma. It is becoming evident that red cells may have a significant influence on
hemostasis and thrombosis; the nature of the effect is related to the flow
conditions.One of the most important
rheological consequences of the multicomponent nature of blood in small vessels
is the migration of red cells towards the core of the flow and commensurate
enhancement of platelet concentration near the boundaries. Several
investigators attempting to simulate the deposition of platelets on biological
and artificial surfaces have found that single-continuum models significantly
under-predict the concentration of platelets near the boundary (see Sorenson et
al. [83, 84], Jordan et al. [85], Anand et al. [86], and Goodman et al.
[87]). This error was attributed to the absence of RBC-induced platelet
margination which enhances the lateral diffusion of platelets towards the wall.
The conclusions of Sorensen et al. and Jordan et al. strongly suggest the
existence of an anisotropic (directionally-dependent) diffusivity of platelets.The mechanism for the above
phenomenon, whereby platelets in flowing blood are propelled towards the
surface, can be understood in terms of their interaction with red blood cells.
Extensive experimental studies dating to the late 1920s (see Fahraeus and
Lindqvist [88, 89]) have demonstrated the phenomenon of lateral migration
of RBCs, resulting in an excess layer of plasma (and platelets) near the wall.
The pioneering work performed by Chien et al. [90–92], Sutera et al.
[93], Blackshear et al. [94, 95], Goldsmith [96–98], and
Hochmuth et al. [99] should be acknowledged. These investigators provided the
earliest insights to understand the microstructural basis for blood rheology.
Most notable is the microscopic “freeze-capture” experiments by
Goldsmith et al. to visualize the margination of red cells in small (65 to
200 mm dia.) cylindrical glass tubes. The associated phenomenon of near-wall
excess of platelets was subsequently studied by Eckstein et al. [100, 101]
and Aarts et al. [102], who empirically revealed the importance of several
independent parameters, including shear rate, hematocrit, the size and
concentration of platelets (or platelet-sized particles), and vessel geometry.In general, blood can be viewed as a suspension and modeled
using the techniques of non-Newtonian fluid mechanics. We, however, assume that
blood is a two-component mixture, composed of the red blood cells (RBCs) suspended
in a (platelet rich) plasma. In the following description, the plasma in the
mixture will be represented byS1 and the RBCs by S2. We
also assume that ϕ represents the concentration of the RBCs, also
commonly referred to as hematocrit. Now, ρ1 and ρ2 are the bulk densities of the mixture
components given by ρ1=γρP,ρ2=ϕρRBC, where ρP is the density of the plasma, ρRBC is the density of the RBCs, ϕ is the volume fraction of the RBCs, and γ is the volume fraction of the plasma. Once the
individual stress tensors are derived (or proposed), a mixture stress tensor
can be defined as Tm=T1+T2,
where T1=(1−ϕ)TPandT2=TR,
so that the mixture stress tensor reduces to that of plasma as ϕ→0 and to that of a RBCs as ϕ→1. Note
that T2 may also be
written as T2=ϕT^R,
where T^R may be thought of as representing the stress
tensor in a reference configuration.We assume that the mixture is saturated, that is,γ=1−ϕ and that the plasma can be represented by the fluid
component of M-R model, where(5.1)p(ρ1)=p(1−ϕ),λ1(ρ1)=λ(1−ϕ),μ1(ρ1)=μ(1−ϕ). And thus (4.1) becomes(5.2)T1=[-p(1-ϕ)+λ(1-ϕ)trD1]1+2μ(1-ϕ)D1, where p is the pressure, μ is the viscosity, and D1 is the symmetric part of the velocity gradient
of the plasma, and λ is the second coefficient of viscosity.Modelling the RBCs is practically more complicated since they
are anisotropic and deformable. The main reason for this difficulty is the
orientation or the alignment of these nonspherical cells. Materials possessing
microstructures, for example, with the internal couples or couple stresses,
were first studied in the early twentieth century by D. Cosserat and F.
Cosserat (Truesdell and Toupin [103]). To study this effect, there are at
least two distinct yet related methods based on continuum mechanics. The first
method is to use an orientation distribution function, whereby one derives
orientation tensors to characterize the behavior of the particles, in this
case, RBCs. The general idea of using orientation tensors to account, in an
averaged sense, for the distribution of fibers in a fluid, in a general
situation, was suggested by Hand [104, 105]. The details of these techniques
are given in Advani and Tucker  [106], Advani [107], and Petrie [108]. The
second method is to use the continuum mechanics theories whereby the
microstructure is in some sense incorporated into the theory, for example, as
is done in the micropolar or director theories (Truesdell and Noll [109]). A
very powerful use of this method is the theory of liquid crystals developed by
Ericksen and later generalized by Hand, Leslie, and others. In this approach, a
unit vector n is introduced as one
of the independent constitutive variables, and as a result the stress tensor
would depend on n and its
derivatives, as well as other important constitutive parameters such as
velocity, velocity gradient, temperature, and so forth, in an appropriate
frame-invariant form.To model the stress tensor for the RBCs, we modify the
constitutive relation derived by Massoudi [110] for an anisotropic granular
media, where the granules were assumed to have a principal direction, denoted
with a unit normal vector n. A general
representation for the stress tensor was given.
Similar constitutive relations have been obtained, for example, by Rajagopal
and Wineman [64] and Rajagopal and Ružička [111] within the context of
continuum mechanics of electrorheological materials:(5.3)T2=a11+a2m⊗m+a3n⊗n+a4(m⊗n+n⊗m)+a5D2+a6D22+a7(m⊗D2m+D2m⊗m)+a8(m⊗D22m+D22m⊗m)+a9(n⊗D2n+D2n⊗n)+a10(n⊗D22n+D22n⊗n)+a11[(m⊗D2n+D2n⊗m)−(n⊗D2m+D2m⊗n)], or(5.4)Tij=a1δij+a2mimj+a3ninj+a4(minj+nimj)+a5Dij+a6Dij2+a7(miDjkmk+Dikmkmj)+a8(miDjk2mk+Dik2mkmj)+a9(niDjknk+Diknknj)+a10(niDjk2nk+Dik2nknj)+a11[(miDjknk+Diknkmj)−(niDjkmk+Dikmknj)], where a1–a11 are scalar functions of the set of  invariants (see
Spencer [112], Wang [113–115], Zheng [116]):(5.5)I1=trD2,I2=trD22,I3=trD23,I4=tr[m⊗m],I5=tr[n⊗n],I6=tr[m⊗n+n⊗m],I7=tr[m⊗D2m],I8=tr[m⊗D222m],I9=tr[n⊗D2n],I10=tr[n⊗D22n],I11=tr[m⊗D2n],I12=tr[m⊗D22n],  where D2=(1/2)[∇v2+(∇v2)T],
and(5.6a)M=m⊗m=gradρ⊗gradρ,(5.6b)Mij=ρ,iρ,j,(5.7a)N=n⊗n,(5.7b)Nij=ninj. This general  equation not only depends on D2 (the symmetric part of the velocity gradient)
and its higher order powers, but also on the density gradient and the RBCs
orientation vector n (a measure of
their anisotropy). There are at least 11 material coefficients in (5.3) which in
some ways have to be specified before a meaningful study can be carried out.
Again, for certain cases, without performing any stability (see Rajagopal et
al. [63]) or thermodynamic analysis, we can gain some information about the
sign of these parameters. Some of the rheological properties can be measured,
for example, using orthogonal rheometers (see Rajagopal et al. [75]). If, for
simplicity, we assume(5.8)a3=a4=a7=a8=a9=a10=a11=0 corresponding to, for example,
a dispersed component such as spherical particles, where there is a degree of
symmetry and the anisotropy of the material does not play a role. However,
density (or volume fraction) gradient is still important. For such a medium, (5.3)
becomes(5.9)T=b11+b2m⊗m+b3D2+b4D22, where b1–b4 are now scalar functions of the appropriate
invariants. Let us furthermore assume(5.10)b1=b1(ρ,trD2,tr(m⊗m)),b2=b2(ρ),b3=b3(ρ),b4=b4(ρ). Now, if  b1 is
given by(5.11)b1=β0(ρ)+β1gradρ⋅gradρ+β2(ρ)trD2, then (5.9) can be written as(5.12)T=[β0(ρ)+β1(ρ)gradρ⋅gradρ+β2(ρ)trD2]1+b2gradρ⊗gradρ+b3D2+b4D22. This equation was derived by Rajagopal and Massoudi [69]. A
special case of this model, with b4=0, has been used extensively by
Massoudi and Rajagopal in a variety of applications (Massoudi et al. [12]). If
we furthermore, for the sake of simplicity, assume that the effects of the volume
fraction gradient are negligible, implying that(5.13)β1=b2=0, then the stress tensor for the RBCs now has the structure(5.14)T2=[βo(ϕ)+β2(ϕ)trD2]1+β3(ϕ)D2+β5(ϕ)D22, where (based on suggestions made by
Massoudi and Rajagopal) the following additional assumptions are made:(5.15)β0=−pϕ,β2(ϕ)=β20(ϕ+ϕ2),β3(ϕ)=β30(ϕ+ϕ2),β5(ϕ)=β50(ϕ+ϕ2). The first  equation, namely, β0=−pϕ,
is different from that suggested by Massoudi and Rajagopal. This equation, in
conjunction with (5.1)a, that is, p(ρ1)=p(1−ϕ),
implies that the total pressure is weighted (distributed) among the two
components according to the volume fraction. This is an accepted assumption in
many two-component theories, and unless there are clear experimental
observations pointing otherwise, this seems to be a reasonable assumption.To account for the shear-thinning effects observed in blood
flow, it is reasonable to assume that the viscosity coefficientβ3 should also depend on Π;
that is(5.16)β3(ϕ,trD2)=β30(ϕ+ϕ2)Πm/2, where(5.17)Π=12tr[2D2]2, when m < 0, the material is
shear-thinning, and if m > 0, it is shear-thickening (see also Pontrelli
[117]). When m = 0, we recover the M-R model. Alternatively, the shear-thinning
effects of blood can also be included in the viscosity term in RBC by assuming
(see Yeleswarapu [118])(5.18)β30=μ∞+(μ0−μ∞)1+ln(1+κγ˙)1+κγ˙, where μ∞ is a term related to asymptotic viscosity
under infinity shear rate, μ0 is a term related to viscosity at zero
shear rate,  γ˙ is the shear rate, and κ is a materials parameter describing the
shear thinning profile. Clearly, other types of viscoelastic fluid models such
as the modified Maxwell or the Oldroyd type fluid models can also be used here.
Substituting (5.15) and (5.16) into (5.14), we obtain, the stress tensor equation for
the RBCs:(5.19)T2=[−pϕ+β20(ϕ+ϕ2)trD2]1+β30(ϕ+ϕ2)Πm/2D2+β50(ϕ+ϕ2)D22. If for simplicity, we define(5.20)λ2=β20(ϕ+ϕ2),μ2=β30(ϕ+ϕ2)Πm/2,δ2=β50(ϕ+ϕ2). Then,(5.21)T2=[−pϕ+λ2trD2]1+μ2D2+δ2D22, which interestingly has the same structure as a Reiner-Rivlin
fluid where the material coefficients are given by (5.20).Finally,
it needs to be mentioned that the RBCs themselves actually form a mixture made
of a central fluid-like region bounded by a viscoelastic solid. Many researchers
have measured the viscoelastic properties of RBCs. However, for the sake of
simplicity, we have decided to treat the RBCs not as a mixture, but as an
isotropic granular-like material.


## 6. A Stress Tensor for Blood

Blood is a complex mixture composed of plasma, red blood
cells (RBCs), white blood cells (WBCs), platelets, and other proteins. A
complete constitutive relation for the stress tensor of the (whole) blood, not
only must capture and describe the rheological characteristics of its different
components, but also must include the biochemistry and the chemical reactions
occurring. To date, no such comprehensive and universal constitutive relation
exists. As mentioned by Anand et al. [86]:
“However, the numerous biochemical reactions that take place leading to the
formation and lysis of clots, and the exact influence of hemodynamic factors in
these reactions are incompletely understood.” In fact, the majority, if not
all, of the papers published on blood characteristics either deal with the
biochemistry of clot formation and other biochemical issues, ignoring
completely the hemodynamic (see, e.g., Kuharsky and Fogelson [119]), or deal
exclusively with hemodynamic or homeostasis and pay no attention to the
biochemical reactions (see Sorensen et al. [83, 84]). Anand and Rajagopal [120]
developed a model for blood that is capable of incorporating platelet
activation. More recently, Rajagopal and colleagues (Anand et al. [86, 121, 122]) have provided a framework whereby some of the biochemical aspects of
blood along with certain rheological (viscoelastic) properties of blood are
included in their formulation.It has been reported that at low shear rates, blood seems to
have a high apparent viscosity (due to RBC aggregation) while at high shear
rates the opposite behavior is observed (due to RBC disaggregation) (see Anand
and Rajagopal [120, 123], Anand et al. [86]). One of the successful models
which has been able to capture the shear-thinning behavior of blood over a wide
range of shear rates is the one proposed by Yeleswarapu [118] and Yeleswarapu
et al. [124]; in this model, the authors proposed a generalization of a three-constant Oldroyd-B fluid. Most of the modeling efforts to use the non-Newtonian
fluid mechanics approach to study blood fall into two categories: (i) those which predict shear-thinning (power law
models or Oldroyd type models), and (ii) models which exhibit yield stress
(Casson model or Herschel-Bulkley type models). However, it is well known that
many suspensions composed of particle (deformable or nondeformable, spherical,
or rod-like) also show normal stress effects giving rise to phenomena such as
die-swell or rod-climbing (see, e.g., Macosko
[125], Larson [126]). Interestingly there is not much experimental work
reported on this subject, that is, whether blood also exhibits normal stress
effects. Massoudi and Phuoc [127] modeled blood as a modified second grade
fluid where the shear viscosity and the normal stress coefficients depend on
the shear rate. It was assumed that blood near the wall behaves as Newtonian
fluid, and at the core it behaves as non-Newtonian fluid.An alternative way is to model blood as we have done here: a
two-component mixture composed of plasma and RBCs. If we are interested in
describing the global behavior of blood, then a stress tensor for blood can be
assumed to be given by(6.1)Tblood=Tplasma(1)+TRBC(2), where by using (5.2) and (5.19), we have(6.2)Tblood=[-p(1-ϕ)+λ(1-ϕ)trD1]1+2μ(1-ϕ)D1+[−pϕ+β20(ϕ+ϕ2)trD2]1+β30(ϕ+ϕ2)Πm/2D2+β50(ϕ+ϕ2)D22. This can be rewritten as(6.3)Tblood=[-p+λ(1-ϕ)trD1+β20(ϕ+ϕ2)trD2]1+2μ(1-ϕ)D1+β30(ϕ+ϕ2)Πm/2D2+β50(ϕ+ϕ2)D22. Now, if we assume that the two velocities are equal, that is,(6.4)v1=v2=v. This refers to an
idealized case where the RBCs are basically being carried away with the
velocity of the plasma. It follows that(6.5)D1=D2=D, and
(6.3) can be rewritten as(6.6)Tblood={-p+[λ(1-ϕ)+β20(ϕ+ϕ2)]trD}1+[2μ(1-ϕ)+β30(ϕ+ϕ2)Πm/2]D+β50(ϕ+ϕ2)D2. If we define(6.7)λblood=[λ(1-ϕ)+β20(ϕ+ϕ2)],(6.8)μblood=[2μ(1-ϕ)+β30(ϕ+ϕ2)Πm/2], then a simple
constitutive relation for blood is obtained(6.9)Tblood=[-p+λbloodtrD]1+μbloodD+β50(ϕ+ϕ2)D2. Now, using a simple viscometer, the apparent blood viscosity
can be measured, where μblood is given by (6.8) indicating that there is a
contribution from the plasma component, and a contribution from the RBCs, appropriately
weighted (see Massoudi [16] for a discussion of this issue. Massoudi [16]
proposed the following hypothesis: “The material properties of one constituent
which appear in the constitutive relation(s) of the other constituent(s),
sometimes referred to as the cross terms, should be multiplied by a function f(ϕ) such that as ϕ→0 or as ϕ→1,f(ϕ)→0 (the simplest case being, ϕ(1−ϕ)), which implies that the cross terms approach
zero, while all other material properties should be weighted according to the
respective volume fraction of that constituent.”). Furthermore,
there is shear-rate dependence from the RBCs and a nonlinear viscosity term,
also a function of the RBCs concentration given by the last term in the above
equation.Obviously, the
form of (6.9) depends both on the form of the constitutive relations for plasma,
(5.2), and for the RBCs, (5.19). Had we assumed an anisotropic representation for
the RBCs, then the stress tensor for blood would have reflected a different fluid-type
characteristic. For example, if in (5.3) we set(6.10)a2=a4=a7=a8=a11=0, which means a flowing anisotropic
material where density gradient does not have an impact on the stress, then we
have(6.11)T2=α11+α3n⊗n+α5D2+α6D22+α9(n⊗D2n+D2n⊗n)+α10(n⊗D22n+D22n⊗n), or(6.12)Tij=α1δij+α3ninj+α5Dij+α6Dij2+α9(niDjknk+Diknknj)+α10(niDjk2nk+Dik2nknj), which is the same as
Leslie-Ericksen (see Leslie [82, Equation (1)], provided that(6.13a)α1=−p, and the α's are functions of(6.13b)nini,Dijninj,DikDkjninj,DijDij,DikDkjDji. Of course, this makes
the problem more complicated in the sense that now, similar to anisotropic liquids,
one needs additional governing (balance) equations and one has to deal with the
unknown additional boundary conditions created by the associated nonlinearities
(see Massoudi [110] for details). The formulation in this section reveals that
the rheological behavior of blood, represented by (6.6), for example, is to a
large extent determined and controlled by the rheological behavior of the RBCs,
as given, for example, by (5.19). We have neglected the deformability of the
cells. In the next section, we will look at a simple boundary value problem.


## 7. Simple Shear Flow between Two Flat Plates

For a simple shear
flow, that is, flow between two horizontal plates  a
distance “h” apart, with the lower plate fixed and the upper plate moving with
a constant speed, the velocity field v and the volume function ϕ are assumed to be of the form (Furthermore, we have assumed that the walls are rigid and nondeformable; however, in reality the vessel walls are flexible and viscoelastic.)(7.1a)v=u(y)i,(7.1b)ϕ=ϕ(y). It then follows that(7.2)D=12(0u′0u′00000),D2=14((u′)2000(u′)20000). Also, notice that(7.3a)trD=0,(7.3b)trD2=12(dudy)2,(7.3c)Πm/2=|dudy|m. Now, using (7.1a), (7.1b), (7.2), (7.3a), (7.3b), and (7.3c) in (6.9), we
find that(7.4a)Txy=12[μblooddudy],(7.4b)Txx-Tyy=0,(7.4c)Tyy-Tzz=[β50(ϕ+ϕ2)](dudy)2. Therefore, we can
see that according to (6.9), blood exhibits only one of the normal stress
differences. If the term β4(ϕ)∇ϕ⊗∇ϕ were kept in the constitutive expression in (5.12),
the model would be capable of exhibiting both of the normal stress differences.Recalling that the
balance of linear momentum for blood as a mixture is given by (3.6), and for the
case of v1=v2=v,
implying D1=D2=D,
we have(7.5)(ρ1+ρ2)DvDt=divTm+(ρ1+ρ2)b, where ρ1=(1−ϕ)ρplasma and ρ2=ϕρRBC,
and(7.6)Tm=Tblood=[-p+λbloodtrD]1+μbloodD+β50(ϕ+ϕ2)D2. For the flowfield assumed by (7.1a), (7.1b), (7.5), with using (7.6), reduces
to the three components in the x, y, and z direction, respectively,(7.7a)12ddy[μblooddudy]+(ρ1+ρ2)bx=0,(7.7b)ddy[−p+β50(ϕ+ϕ2)4(dudy)2]+(ρ1+ρ2)by=0,(7.7c)(ρ1+ρ2)bz=0, where bx,by,andbz are the components of the external body force,
and μblood is given by (6.8). Thus, we see that a motion
of the form of (7.1a) and (7.1b) is only possible if the z-component of the body force
field is zero.Equations (7.7a) and (7.7b) form a system of two coupled
second-order nonlinear ordinary differential equations, and in general have to
be solved numerically. Appropriate boundary conditions are needed in order to
have a well-posed problem. For the simple shearing motion assumed by (7.1a) and (7.1b), we
assume(7.8a)u(0)=0,(7.8b)u(h)=V, where “h” is the distance between the two plates. It is
possible that the no-slip condition may not be appropriate for all cases, as
the RBCs might roll or slip at the boundaries. And for this special case, where
the gradients of volume fraction are ignored, we only need one boundary
condition for ϕ:(7.9)ϕ(0)=ϕ0, where ϕ0 is a constant (for different types of boundary
conditions, see Langtangen and Munthe [128], and Massoudi [41]).To obtain an analytical solution (closed form solution) to
the above (7.7a), (7.7b) under some idealized conditions,  we furthermore assume(7.10)bx=by=bz=0. Then, (7.7a) can be integrated once to  give(7.11)[μblooddudy]=C1, where C1 is a constant. Recalling that (7.3c), and (6.8), we have(7.12)μblood=[2μ(1-ϕ)+β30(ϕ+ϕ2)Πm/2]=2μ(1-ϕ)+β30(ϕ+ϕ2)|dudy|m. Substituting this in (7.11), we have(7.13)[2μ(1-ϕ)+β30(ϕ+ϕ2)|dudy|m]dudy=C1. Equation (7.13) can be integrated once to give us(7.14)−p+β50(ϕ+ϕ2)4(dudy)2=C2, where C2 is a constant. Now, it can be shown that the
above equations admit a solution of the form(7.15)u=αy,ϕ=k, where α and k are constants. However, the system of (7.7a)
and (7.7b) is nonlinear and might admit additional solutions.


## 8. Comments

In this paper, we have discussed, based on the classical mixture
theory and the approach taken by Massoudi and Rajagopal, a framework for modeling
the rheological behavior of blood. The proposed constitutive relation depends
on the form of the stress tensors for the plasma and the red-blood cells. In general,
the RBCs are assumed to behave as an anisotropic density-gradient-dependent viscous
fluid. As such the equations are highly nonlinear. After making many
simplifying assumptions, a relatively simple constitutive relation for the
stress tensor for blood is obtained (see (7.6)); its form is very similar to a
Reiner-Rivlin fluid, where the shear viscosity coefficient is not only a
function of the shear rate, but also of the concentration. It
is noted that we have only discussed the development of a model and that
specific boundary value problems need to be solved to test the efficacy of the
model. A simple shear flow is studied, and an exact solution is obtained for a
very special case; for more general cases, it is necessary to solve the
nonlinear coupled equations numerically. Furthermore, it should be
mentioned that for higher order or higher gradient theories, assigning boundary
conditions for certain terms, which appear in the governing equations, is a
difficult task. Quite often these boundary conditions are not derived from
first principles; instead they are given as ad hoc assumptions, or they are
simply specified as mathematical conveniences. Sometimes experiments have been
used successfully to specify these necessary additional boundary conditions.Furthermore, it goes without saying that the model
developed here is only appropriate for a healthy human, and it does not capture
any blood disorder.  It has also been shown that with regular exercise and
physical training certain characteristics of blood can change, and blood
undergoes what is known as “fluidification” (see Ernst and Matrai [129],
Wannamethee et al. [130]).To include
the formation and growth of clots, and lysis of blood cells in blood, in
general, the reaction-convection-diffusion equations are to be solved in
conjunction with the balance laws for mass, linear and angular momentum, and
energy (for each component). Although we have ignored the biochemical effects
of blood in this paper, in principle, the theory is amenable to  extension (see
Anand et al. [86]).



---

*Source: 579172-2008-09-08.xml*