THREEPHASE HYBRID MIXTURE THEORY FOR SWELLING DRUG
DELIVERY SYSTEMS
by
Tessa F. Weinstein
M.S., University of Colorado at Denver, 2002
B.S., Metropolitan State College of Denver, 1999
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
2005
This thesis for the Doctor of Philosophy
degree by
Tessa F. Weinstein
has been approved
by
Tom Russell
Ron Rorrer
Weinstein, Tessa F. (Ph.D., Applied Mathematics)
ThreePhase Hybrid Mixture Theory for Swelling Drug Delivery Systems
Thesis directed by Associate Professor Lynn S. Bennenthum
ABSTRACT
Drug delivery systems are composed of a polymeric' carrier, such as hydrox
ypropyl methylcellulose (HPMC), and a drug or active agent, such as naproxen,
the drug found in brands names such as Aleve. When such a device is ingested,
the polymer swells and the porespace increases in size, the drug diffuses out
of the device. Some of these devices are swellingcontrolled, meaning that the
rate of swelling is the controlling factor in release kinetics. Modeling swelling
controlled polymeric drug delivery systems presents unique problems. As the
solid drug undergoes phase transfer, the density of the solid phase, composed
of polymer and drug, changes. As such, the usual incompressibility condition,
that the density of the solid phase does not change in time, no longer applies.
To overcome this difficulty we modify a twophase mixture theoretic model to
a threephase model where the polymer and drug are modeled as separate solid
phases, and the liquid phase remains water and drug. A twoscale and three
scale theory for such a threephase model is developed. First, it is shown how a
new choice of independent variables yields a physically meaningful interpretation
of the solid phase stress, pressure, and Terzaghi stress tensors. We present the
threephase theory modeling the polymer as viscoelastic, the drug as elastic, and
the liquid as a viscous fluid. We present constitutive theory, nonequilibrium re
sults, equilibrium results, nearequilibrium results, as well as Darcys and Fick's
laws for both the twoscale and threescale scenarios. Appropriate simplifica
tions that are specifically applicable to swellingcontrolled drug delivery systems
are presented and used to combine the Darcy's and Ficks laws with the con
servation of mass equations to obtain a system of equations which models the
transport of the drug in terms of liquid volume fraction and the concentration
of the drug in the liquid phase.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
ACKNOWLEDGMENT
First, and foremost, I would like to thank Lynn S. Bennethum, my advisor, for
her neverending patience and direction. I only hope that I can guide my first
student as well as she has guided me through this. Id also like to thank John
H. Cushman for the opportunity that he provided me during my last year as a
Ph.D. student, which was instrumental to me finishing on time.
CONTENTS
Figures ............................................................ xi
Tables............................................................. xii
Chapter
1. Introduction..................................................... 1
1.1 Terminology...................................................... 1
1.2 Previous Work ................................................... 6
1.3 Thesis Outline .................................................. 9
2. Overview of Hybrid Mixture Theory............................... 12
2.1 Microscale Governing Equations.................................. 13
2.2 Averaging Procedure............................................. 14
2.3 Mesoscale Field Equations....................................... 17
2.3.1 Conservation of Mass......................................... 17
2.3.2 Linear Momentum Balance ..................................... 20
2.3.3 General Case................................................. 22
2.3.4 Angular Momentum Balance .................................... 23
2.3.5 Conservation of Energy ...................................... 28
2.3.6 Entropy Balance.............................................. 31
2.4 Standard Assumptions and Resulting Entropy Inequality........... 33
2.5 Choosing Constitutive Independent Variables and Exploiting the En
tropy Inequality................................................ 35
vii
2.5.1 Example: Elastic Solid......................................... 36
2.5.2 Example: Thermoviscoelastic Fluid .............................. 39
3. A New Choice of Independent Variables............................. 44
3.1 Motivation....................................................... 44
3.2 Constitution..................................................... 48
3.3 The Entropy Inequality .......................................... 50
3.3.1 Two Phase Pressures.............................................. 56
3.4 Novel Results.................................................... 61
3.4.1 Solid Phase Stress .............................................. 61
3.4.2 Darcys Law...................................................... 63
3.5 Transport Models................................................. 67
3.5.1 Transport Equation in Terms of Pressures......................... 68
3.5.2 Transport Model in Terms of the Chemical Potential............... 72
3.5.3 Model Comparison................................................. 74
3.6 Conclusions...................................................... 75
4. ThreePhase HMT for Swelling Porous Systems ........................ 76
4.1 Motivation....................................................... 76
4.2 Constitution..................................................... 78
4.3 Entropy Inequality............................................... 80
4.3.1 Three Phase Pressures............................................ 84
4.3.2 NonEquilibrium Results.......................................... 87
4.3.3 Equilibrium Results ............................................. 90
4.3.4 Removing Nth Component Dependencies...................... 95
4.3.5 NearEquilibrium Results......................................... 99
viii
4.3.6 Darcys Law..................................................... 103
4.3.7 Ficks Law ..................................................... 104
5. TwoScale Polymeric Drug Delivery Systems........................ 105
5.1 Bulk Transport.................................................. 105
5.1.1 Species Transport............................................... 109
5.2 Discussion...................................................... Ill
6. ThreeScale HMT for Swelling Porous Systems...................... 113
6.1 Macroscale Notation and Averaging Procedure..................... 113
6.2 Macroscale Field Equations...................................... 117
6.2.1 Conservation of Mass............................................ 117
6.2.2 Conservation of Linear Momentum................................. 119
6.2.3 Conservation of Energy ......................................... 120
6.2.4 Entropy Production ............................................. 123
6.3 Standard Assumptions and Resulting Entropy Inequality........... 124
6.4 Constitution.................................................... 126
6.5 Entropy Inequality.............................................. 129
6.5.1 Macroscale Pressures............................................ 134
6.5.2 NonEquilibrium Results......................................... 137
6.5.3 Equilibrium Results ............................................ 140
6.5.4 Removing Nth Component Dependencies............................. 147
6.5.5 NearEquilibrium Relations...................................... 151
6.5.6 Darcys Laws.................................................... 155
6.5.7 Ficks Law ..................................................... 156
7. ThreeScale Polymeric Drug Delivery Systems...................... 158
IX
7.1 Vicinal Fluid .................................................... 159
7.2 Bulk Fluid, B..................................................... 161
7.3 Fluid Phase, C ................................................... 162
7.4 Bulk Transport.................................................... 163
7.5 Species Transport................................................. 168
7.6 Discussion........................................................ 171
8. Discussion and Avenues for Future Work.............................. 173
Appendix
A. Mesoscale Appendix ................................................. 175
A.l Nomenclature...................................................... 175
A.2 Bulk Phase Definitions............................................ 181
A. 3 Identities Needed to Obtain Entropy Inequality 2.82.............. 183
B. Macroscale Appendix................................................ 184
B.l Nomenclature...................................................... 184
B.2 Bulk Phase Definitions............................................ 189
References............................................................. 192
x
FIGURES
Figure
1.1 Scales of Observation............................................... 2
1.2 Different types of polymer molecules................................ 3
1.3 Observed Fronts .................................................... 6
6.1 Macroscale REV.................................................... 114
xi
TABLES
Table
2.1 Quantities for Microscopic Field Equations
1. Introduction
Porous materials consisting of a porous solid with fluid filled or gas filled
pores appear in a wide variety of applications. In particular, polymers have
become increasingly important in technological industries. Polymers have a
myriad of applications including but not limited to: construction, agriculture,
drug delivery, and food stuffs. The focus of this thesis will be on the swelling
polymers polymers that swell (shrink) upon wetting (drying) used in drug
delivery applications.
Swelling polymers exhibit a hierarchy of scales due to their complex porous
structure. Herein, we develop both a twoscale and a threescale theory for
swelling polymer systems. We call these scales the microscale, mesoscale, and
macroscale. Figure 1 depicts the relevant scales of polymer particle immersed
in a bulk fluid.
1.1 Terminology
Before we discuss the complex interplay of forces at work in these systems,
a basic understanding of polymer terminology, characteristics, and behavior is
needed. Although there is still some debate on the subject, polymers generally
have molecular weights of 25,000 g/mol or larger, and increased chain length
means increased entanglement of polymer chains. Monomers are the building
blocks of polymers, e.g. amino acids or sugars. Monomers link together with
other molecules of the same or different type to form polymers. A homopolymer
is made up of a single monomer, whereas a copolymer is made up of two or more
1
Macroscale
bulk fluid
Mesoscale REV
mesoscale polymer
polymer
bulk fluid
polymer
vicinal water
Mesoscale
Microscale
Figure 1.1: Scales of Observation
2
Hr.iiKhed
Network
Figure 1.2: Different types of polymer molecules
types of monomers. A repeating unit is a segment of a macromolecule that forms
the basic unit of the macromolecule (excluding the ends). In other words, we
can form a complete polymer by linking an adequate number of repeating units
together. A linear polymer is one in which every repeating unit is attached to
exactly two others. A branched polymer is one in which repeating units are not
linked exclusively in linear fashion. See figure 1.1 for a schematic representation.
A network polymer is formed by chemically linking together linear or branch
polymers; this process is referred to as crosslinking. This same process is known
as vulcanization for rubbers. See [45, 55, 58, 54].
Polymers can be broadly classified into two groups: amorphous polymers,
which can be further classified as crosslinked or uncrosslinked. and crystalline
polymers. Amorphous portions of a polymer are the areas of the polymer with
long chain length where the polymer tends to coil about itself randomly with no
apparent order. Unvulcanized natural rubber is an example of an uncrosslinked
amorphous polymer. Vulcanization of natural rubber introduces covalent sulfur
bonds so that the chains are chemically bonded to one another. Crystalline
polymers exhibit a crystalline structure, however, most of them contain regions
of amorphous material.
3
Polymers exhibit a wide spectrum of behavior dependent on various proper
ties such as degree of crystallinity, crystal size, crosslinking, molecular weight of
the polymer, and surrounding temperature, all of which can affect diffusive be
O
havior. In systems with fluid filled pores on the order of 200 500 A, molecular
diffusion occurs due to concentration gradients. In less porous systems molecu
lar diffusion of the drug may happen throughout the entire polymer phase, and
dissolution of the polymer matrix may be necessary to accomplish release of the
active agent (drug). Some polymers are so highly sensitive to temperature that
certain properties change almost discontinuously at a temperature, Tg, called the
glass transition temperature, [18]. Such polymers are called thermosets. Below
the glass transition temperature polymers are solid, hard, even brittle and are
considered to be in the glassy state. Amorphous polymers in the glassy state are
sometimes called amorphous liquids or supercooled liquids. It is well known that
amorphous polymers below Tg are not in thermodynamic equilibrium; they still
flow, but the time scale for observing creep and flow is very long, as is the case
with common glass. Above Tg polymers enter the glassrubber transition, where
the polymer softens. As the temperature increases polymers will enter a rub
bery plateau, rubbery flow state, and finally a viscous flow state [54]. Diffusion
properties of a polymer change drastically as the polymer nears Tg [36],
Hydroxypropyl methylcellulose (HMPC) is the prevalent polymer carrier
used in oral drug delivery systems, [47]. Because of the hydrophilic nature of
HMPC, it displays strongly timedependent behavior when exposed to water,
swelling considerably, and eventually dissolving. In a recent experiment pre
formed by Colombo et ah, [20], cylindrical disks containing HMPC and varying
4
initial loadings of flurorescein sodium (a water soluble dye used as a detection
agent in visible images) were clamped between transparent Plexiglass slabs, in
troduced to water filled vessels, removed at different times, and photographed.
The results showed that, depending on the initial loading of the active agent,
three distinct fronts were observed. Consider figure 1.1. The swelling front is
denoted by S. As water penetrates the polymer matrix it acts as a placticizer
by increasing the free volume, reducing the glass transition temperature of the
system, Tg. Eventually, Tg will equal the experimental temperature, T, the
polymer enters the glass transition, as discussed above, and swells. Thus, inside
front S, Tg > T and the polymer remains in a glassy state, and outside front
S, Tg < T and the polymer is in a swollen, rubbery state. As water penetrates
the matrix, drug dissolves and diffuses out of the material. The rate at which
the drug diffuses out of the material is directly proportional to the amount of
water that has penetrated the material. At high enough drug loadings a dif
fusion front, D, is observed because when the maximum solubility of the drug
is reached both dissolved drug, and nondissolved drug (which is not available
for diffusion) exist within the polymer matrix. Lastly, an erosion front, E, is
observed due to polymer disentanglement caused by the molecularlevel snake
like motion of the polymer chains (reptation) [47]. The effects of the erosion
front on release kinetics depends largely on whether or not all of the drug has
diffused out of the polymer matrix by the time the erosion starts to occur. The
effects of the erosion front become of greater concern in systems where the drug
has a low solubility. Then the dissolution of the polymer matrix becomes the
controlling factor in release kinetics.
5
Figure 1.3: Observed Fronts
1.2 Previous Work
Traditional models of drug delivery systems begin with a concentration form
of Ficks law of diffusion, such as
(1.1)
where q is the mass flux and C is the concentration of the active agent. De
pending on the type of system modeled and geometry, Ficks law is modified
accordingly. Classical models use the effective diffusion coefficient
Deff = Diw, (1.2)
where DIW is the diffusion coefficient of the drug when the pores are filled with
a solution of the drug, e is the porosity, and r is the tortuosity. The diffusion
6
coefficient can be modified by considering it as the mutual diffusion of the drug
in the polymer, which measures the change in the concentration of a species
from its average concentration with respect to time, the diffusion of the drug
alone, or making it concentration dependent. The primary problem with these
variations is that they are of limited applicability since they only apply to Fickian
diffusion phenomena [43]. For swellingcontrolled drug delivery systems, it is well
known that nonFickian diffusion occurs when polymers enter the glassrubber
transition [36], Various authors have expanded the Fickian model to account
for anomalous behavior, but these extensions are largely heuristic [43]. More
rigorous mechanistic models exist, such as that proposed in [48]. However, they
assume a priori that the diffusion coefficient of the drug is a function of the
water content (liquid volume fraction) of the system to capture the dependency
of the rate of drug dissolution on water content. Our analysis will show that a
gradient in volume fraction should be used to account for this behavior.
The author notes here that it is common in the drug delivery literature for
authors to refer to equations for the fluid phase penetration into the polymer as
Ficks law. The bulk fluid is made up of species including: water, active agent,
and perhaps dissolved polymer. We will adopt the hydrology convention and
use Ficks law to refer to diffusion of these species within the fluid phase, and
use Darcys law to refer to penetration of the bulk fluid into the polymer matrix.
A variation on the traditional approach is to treat coefficients of Ficks
and Darcys law as stochastic processes using the KarhunenLoeve expansion,
[27], However, this increases the dimension of the problem by treating the
random aspect of the problem as a new dimension and it is not entirely clear
7
how the resulting coefficients relate back to the physical problem. Lustig et
al. [36], address the problem using continuum thermodynamics. However, they
do not have a variable that directly accounts for the moisture content of the
polymer. Low [34] found that the swelling pressure in smectitic clays (specifically
montmorillonite), is highly dependent on moisture content. We find it reasonable
to assume that the same will be true of the swelling of polymers.
Hybrid mixture theory (HMT) has recently been applied to polymeric and
biopolymeric systems by Singh et al. [49, 50, 51]. Hybrid mixture theory has
all of the advantages of classical mixture theory plus the added advantage that
mesoscale variables can be directly related to microscale counterparts. This last
method is the context in which our model is systematically developed.
HMT consists of averaging the microscopic conservation and balance equa
tions (mass, momentum, energy) to obtain macroscale analogues, then apply
ing Bowens continuum theory of mixtures to obtain constitutive relationships
[15, 16, 17] by exploiting the entropy in the sense of Coleman and Noll [19].
The approach was pioneered by Hassanizadeh and Gray in a series of papers
[29, 30, 31], From 1994 to present HMT has been successfully employed to
model swelling and shrinking behavior in gels, food stuffs, and colloidal systems
where phase interactions play an important role in the mesoscopic and macro
scopic behavior [2, 7, 8, 12, 39, 40, 41, 42, 49, 50, 51]. Most of these works
assumed that the solid phase is composed of an elastic material and a viscous
liquid phase at the microscale. Using HMT the resulting macroscopic models
exhibited viscoelastic behavior due to phase interactions. In [50, 51], a model
was developed in which the solid itself is modeled as viscoelastic, and the result
ing bulk transport equation is an integrodifferential equation in terms of the
volume fraction which has been successfully used to model the drying of food
stuffs.
Some the aforementioned papers [2, 7, 8] contained the incorrect result that,
at equilibrium, the chemical potential of two phases are equal to each other. This
contradicts the classical thermodynamic result that, at equilibrium, the chemical
potential of a single constituent in two phases is equal. In 1994, Achanta et al.
correctly employed HMT with the additional axiom of equipresence [1], which
requires that before exploiting the entropy inequality it is assumed that each
phase is composed of the same N constituents. Only after the entropy inequality
has been exploited can the concentration of certain constituents be set to zero.
In [1], the authors derive the macroscale field equations for each phase and
interface of a threephase, multiconstituent media.
The purpose of this thesis is many fold. One of the goals of this thesis is
to reexamine the assumptions used in [49] and in the process develop ther
modynamic quantities which are physically meaningful. However, the primary
aim of this thesis is to extend the results of Singh et al. [49, 50, 51] to in
clude constituent transport so that we can accurately model swellingcontrolled
drug delivery systems using a polymer carrier such as HMPC. Doing so presents
some unique problems which are outlined below and discussed in great detail
throughout this thesis.
1.3 Thesis Outline
For the sake of completeness, Chapter 2 will give an overview of HMT. We
will start by summarizing the averaging procedure and giving a detailed deriva
9
tion of the mesoscale conservation and balance equations. We will discuss some
standard assumptions, give the resulting entropy inequality, discuss how the in
dependent variables are chosen, and clarify the process of exploiting the entropy
inequality in the sense of Coleman and Noll [19] by using two demonstrative
examples.
A common result of using HMT is that the thermodynamic definition of solid
phase pressure does not coincide with the actual physical stress exerted on the
solid phase. In Chapter 3 we present a solution this problem. First, we present
the constitutive assumptions for a twoscale, twophase system consisting of a
viscoelastic solid phase and a viscous liquid phase, motivate our new choice of
independent variables, and derive the entropy inequality using our new choice
of independent variables. We will present only the novel results which give rise
to a physically meaningful interpretation of the solid phase stress, pressure, and
Terzaghi stress tensors. To show the usefulness of this new set of independent
variables, we will derive a transport equation analogous to that found in [51].
Using Darcys law and FloryHuggins theory we will develop another transport
equation specifically applicable to polymers.
A fewT problems arise when trying to apply twophase HMT to drugdelivery
applications. In [50], Singh et al. assume that the solid phase is incompressible,
which is perfectly reasonable for the system under consideration in that pa
per since it was aimed at systems in which no mass transfer takes place between
phases. However, when considering mass transfer of constituents between phases
the incompressibility condition for the solid phase no longer holds since the den
sity of the phase changes with time as drug leaves the solid phase. Furthermore,
10
experimental literature often refers to the volume fraction of drug present in the
system, and unfortunately twophase HMT has no way straight forward way of
expressing this variable. Using three phases instead of two phases has the ability
to overcome both of these problems. In Chapter 4 we will present the constitu
tive theory for a twoscale, threephase system consisting of a viscoelastic solid
(polymer), an elastic solid (drug), and a viscous fluid (bulk fluid), keeping the
theory as general as possible. We will then derive the entropy inequality for this
system, and present nonequilibrium, equilibrium, nearequilibrium results, as
well as Darcys and Ficks laws.
In Chapter 5 we will present the model assumptions specifically applicable
to swellingcontrolled drug delivery systems, using them to combine the consti
tutive equations with the mesoscale conservation and balance equations from the
previous chapter to derive a novel bulk fluid transport equation and constituent
transport equation. Additionally, we will discuss the experimental data needed
to solve these equations.
In Chapter 6 we will present the constitutive theory for a threescale, three
phase system consisting of a particle composed of viscoelastic solid, an elastic
solid, a viscous fluid and two additional bulk fluids. We will then derive the
entropy inequality for this system, and present nonequilibrium, equilibrium,
near equilibrium results, as well as Darcys and Ficks laws. Chapter 7 will
mirror Chapter 5 but for the threescale system. We will discuss avenues of
further research in the final chapter.
11
2. Overview of Hybrid Mixture Theory
The purpose of this chapter is to give a comprehensive overview of the
averaging procedure used in HMT for any twoscale, multiconstituent, multi
phase material. The two scales are herein called the microscale and mesoscale.
At the microscale one can distinguish between phases. It is at this scale that the
field equations (conservation of mass, momentum balance, etc.) are known to
hold and properties such as density, velocity, and mass are clearly defined. The
mesoscale is order of magnitudes larger than the microscale; at the mesoscale
one can no longer distinguish between individual phases. Because an averaging
process is performed to obtain an analogue on this larger scale, properties such
as density and velocity are now viewed as somewhat blurred. Herein we develop
one such technique.
We assume that the material we are modeling has negligible interfacial ef
fects; that is, the interface has no thermodynamic properties and is massless.
As such, no constituent present gains or loses mass, momentum, or energy when
crossing an interface. This will place special restrictions on each of the field
equations and will be discussed in further detail in subsequent sections. Second,
the material we are modeling has a representative elementary volume (REV);
that is, a volume for which averaged properties will remain the same if the REV
is made slightly larger or smaller. In addition we require that the REV size and
shape remain the same for all space and time. Such an REV does not exist if the
material under question is too heterogeneous. For a more detailed discussion of
12
the existence of such an REV see [5, 21]. It is important to remember that for
the current discussion, the following theory is applicable to any material meeting
these requirements.
2.1 Microscale Governing Equations
Each phase is denoted by small Greek letters (a, /?, 7), and species (or
constituents) are denoted by j, j = 1,N. In general, a subscript Greek
letter indicates a macroscale quantity from that phase. Superscript minuscules
indicate the constituent, so that, e.g. vaj is the macroscopic velocity of con
stituent j in the aphase. A carrot over the symbol, is used to emphasize
that the quantity represents a transfer from either another phase or from other
constituents. A complete nomenclature for the microscale to mesoscale averag
ing and mesoscale variables can be found in Appendix A. For the purposes of
simplicity and brevity we assume that all interfaces are massless and have no
thermodynamic properties. Interfacial effects can be included by following the
approach provided by [28].
Additionally, we assume that there are no internal surface discontinuities,
meaning that each phase is the union of several isolated simply connected vol
umes, [29]. Using the notation of Eringen, [24], the constituent, microscopic
field equations for a given phase, a, can be stated as
+ V {p>vjp3) V i3 (jp = p>G3 + p>p (2.1)
where pi is the mass density, ip3 is the mass average (over the phase) thermody
namic property of constituent j, vJ is the mass average velocity vector, V is the
flux vector, f3 is the body source, G3 is the net production, and ip accounts for
13
Table 2.1: Quantities for Microscopic Field Equations
Quantity ijj i f ip G
Mass 1 0 0 r 0
Linear Momentum V t 9 i + rv 0
Angular Momentum r x v r x t r x g r x [i + rv) rri 0
Energy e+lv2 t v + q g v + h Q+i vJr?{e + \v2) 0
Entropy V 4> b rj + rv A
the influx of ip from all other constituents (e.g. due to chemical reactions). If the
medium consists of only one constituent then ip is equal to zero. This equation
holds on the microscale for mass, linear and angular momentum, energy, and
entropy. Table 2.1 lists the quantities used for each field equation.
In Table 2.1, t is the secondorder stress tensor, g is the external supply of
momentum (gravity), r is the position vector referenced to a fixed coordinate
system, e is the internal energy density function, h is the external supply of
energy, q is the heat flux, rj is the entropy density, is the entropy flux, b is the
external supply of entropy, and A is the entropy production.
2.2 Averaging Procedure
The averaging procedure used in HMT is based on works of various au
thors and was developed at approximately the same time [3, 38, 56, 57]. While
many different methods are available [21], we choose to use the one which is
computationally the simplest where field equations are averaged via weighted
integration. Here we use the indicator function of the aphase as the weight,
14
and treat the averaged quantities resulting from the weighted integration as dis
tributions. This allows us to bypass the difficulties of defining the derivative of
averaged quantities that result from the weighted integration [46, 44]. Addition
ally, the weighting function used here may result in averaged quantities that do
not correspond to physical quantities measured. This problem can be overcome
by choosing a weighting function that represents the experimental apparatus
used to measure physical properties [21].
Let 5V denote the REV, 6Va denote the portion of the ophase within 8V,
and 5Aag denote the portion of the a(3 interface within 5V. It is assumed that
8Va and 5Aag are isolated simply connected regions. Expressing the magnitude
of 8V by <5V, the volume fraction can be written as
E(xi) = M (2'2)
so that
X>Q = ! (23)
Q
Letting r and x denote the position vector and the centroid of the REV, respec
tively, r can be written
r = x + Â£, (2.4)
where Â£ is the local coordinate referenced to the centroid of the REV and varies
over all of 8V. The indicator function for the aphase is given by
{1 if r e 8Va
0 if r e SVp,
Then
\SVa\(x,t)= [ ja(x + Â£,Â£) dv(Â£) (2.5)
Jsv
15
represents the magnitude of the volume 6V in the ophase. Following Has
sanizadeh and Gray [29], we make the following definitions:
/ P*(r,t)'ya{r,t)dv(Â£) (2.6)
\oVQ\ Jsv
is the average mass over 5VQ,
(^)(x,t) = ^ / ip3{r,t)fQ{r,t) dv{Â£) (2.7)
IoFqI jsv
is the volume average property of ip3, and
FJQ(a:,t) = tt1 [ g,'(r!t)^(r,t)7a(r,i) dv(Â£) (2.8)
p> \6Va\ Jsv
is the mass average property of yP. We would like the mesoscale field equations
to be analogues of the microscale field equations. To ensure this we will apply
the following theorem that allows us to interchange the order of differentiation
and integration. This result is due to Cushman [22],
Theorem 2.1 If wa>3 is the microscopic velocity of the interface a/3 and na
is the outward unit normal vector of dV0 indicating the integrand should be
evaluated in the limit as the a(3interface is approached from the aside then
w\Ldiiiv{i)=!i[w\LMv(i)
 E Twif fws "v( <29>
1 1JA
w\Lv,idvi=v[\kLMi).
(2.10)
We now have everything we need to average equation (2.1) from the mi
croscale to the mesoscale.
16
2.3 Mesoscale Field Equations
In this section we will derive the mesoscale field equations. We begin with
the conservation of mass, linear momentum balance, angular momentum bal
ance, energy balance, and lastly entropy. Corresponding bulk equations are
derived for each field equation and restrictions that result from the assumption
that the interface has negligible thermodynamic properties are given.
2.3.1 Conservation of Mass
Substituting the appropriate quantities from Table 2.1 into equation (2.1)
the conservation of mass at the microscale is
(pj) + V (p3v3) = p3?3. (2.11)
Formally, we multiply by the indicator function, integrate over the mesoscopic
REV, and divide by the magnitude of the REV, 5V. Using Theorem 2.1,
definitions (2.6)(2.8), and substituting back into equation (2.11), we have
d / nra\ __ / :arc
Qt (^ V ) + V {Â£ p> UJ
1
Y j ^(w010! vJ) nada + Â£api*?3 . (2.12)
Qy
Now, using ^2., the material time derivative, which is given by
DQ] d n _
Dt dt+
(2.13)
we obtain the mesoscopic mass balance for constituent j in phase a
Da> {eapai)
Dt
+ eQpa3{V va>) = +?*3
(2.14)
17
In the equations above, vj is the microscopic velocity of constituent j, wa@j is
the velocity of the jth constituent in the a,3 interface, and na is the outward
unit normal vector of 5Va. Complete nomenclature is given in Appendix A. Our
motivation in defining the mesoscopic variables is twofold. First, we would like
the mesoscopic variables to coincide with actual physical variables that can be
measured via mesoscale experiments. Second, we would like the definition of the
mesoscopic variables to be as consistent as possible with their microscale coun
terparts. Unless otherwise stated, the mesoscopic variables have the following
definitions:
Pa = f*
(2.15)
is the average mass over 6Va,
v j = vJ
(2.16)
is the mass averaged velocity,
T03 =  v>) na da (2.17)
is the net rate of mass gained by constituent j in phase a from phase /3, and
r3' = eapa>?J
(2.18)
is the rate of mass gain due to interaction with other species within the same
phase.
Bulk phase variable definitions are not always intuitive and are defined so
as to preserve the form and interpretation of the mesoscale equations. To obtain
18
the bulk phase counterpart for the conservation of mass we make the following
definitions:
N
Pa = Y,paj (219)
3=1
is the mass density of phase a, and
Cai =
paj
Pa
(2.20)
is the concentration of the jth constituent in the aphase. All other definitions
of bulk phase variables are given in Appendix A. Using equations (2.19) and
(2.20) in equation (2.12) we have:
^{eapaCa>) + V (eapaCava) = +
(5
(2.21)
Summing over j = 1,... ,N we obtain:
(Â£Y) + V(eV*a) = (222)
Defining ^ as + va V, we obtain the bulk phase counterpart for the con
servation of mass,
Da(sapa) W(V. where the following restrictions apply: (3^a (2.23)
N
E^0 Va and (2.24)
E ri= a,i3^oc (2.25)
The meaning of equation (2.24) is that the gain of mass of the bulk phase due
to chemical reactions alone must be zero. Equation (2.25) says that the rate of
19
mass gained by phase a from phase 0 is equal to the rate of mass gained by
phase 0 from phase a, i.e. no mass is lost in the interface.
2.3.2 Linear Momentum Balance
Substituting the appropriate quantities from Table 2.1 into equation (2.1)
the conservation of linear momentum at the microscale is
(pv') + V (p^vjvJ) V t7 fPgJ = p?% + prv^. (2.26)
Using the same method as for the conservation of mass equation we obtain:
J^eV'v0') V (Â£a({tj)a + paiva>va> p^WQ))
+ V (EapajVa:>Vaj) Â£apaigaj
= E IaTTT / ^ + piv3(w0i v3)} n da
8^a ' '
_________
+ eapa*(i + ?V), (2.27)
where g"1 = g1 is the mass average external supply of momentum and all other
variables retain the meaning they were given in the previous section. Addition
ally, we define the mesoscopic stress tensor, also known as the Cauchy stress
tensor, as:
ta> = (tJ)a + pava>va> pa^vH0. (2.28)
Now, subtract vQj times equation (2.12) from equation (2.27) to obtain
r\
Â£apa>^t(vai) + Â£apava>V (d>) V (Â£ata>) Â£apajga
T7777 [ W + p1v\walii vj)] na da
r M Ja
f pj(wat3j vj) na da
J A
+ Â£apa> (?   v. (2.29)
E
VaJ
W\
20
Using (2.13) we obtain the mesoscopic linear momentum balance for constituent
j in the aphase:
Â£apaJ
~Dt
("') V {enV>) Â£apn>g] = T3 + J2 T3,
(2.30)
where we have made the following definitions:
T =
3
\svc
 v3
[tj + p?vi(wa!3i vj)\ na da
[ Pi'i
J A
wa0i v3) n da
(2.31)
represents the effect constituent j of phase /3 has on the rate of change of me
chanical momentum of the same constituent in phase a, and
i 1 =Â£apaj(i +rJib )
(2.32)
is an exchange term that takes into account all gain of momenta due to the
presence of other species but not due to chemical reactions.
Again, we would like to obtain the bulk phase counterpart. Beginning with
equation (2.27) and using (2.20) we have
d
{Â£apaC3va3) V (e(t> Pa>va3'a vta)) sapaCajgaj
+2V{Â£Qpa3C3va3 va)~ V{Â£apaCaiva va)
= f [tJ + p^vj(wai}3 uJ)] na da
+Â£apa3{? +Â¥&). (2.33)
21
Subtracting va times equation (2.21) we get
r\ r\
^{sQpaCa>'v>) va~(sapaC^) (2.34)
V (Â£a(taj paiVaj'a Vai,a))
eapaCaigai + V(2eapajCajvaj va)
V{eapaCa>va va) vaV (Â£paca>vai)
= Y Ta0j + Y^vaj + + ?3vaj ~ 53 (235)
3^cl
Summing from j = 1,..., TV we obtain
Dava Â£ap n V (Â£ata) H Dt v 1 eaPaga = Yfl (2.36)
where the following restrictions apply:
N +raivi,a) = 0 Va,and (2.37)
Â£ (r? + '8?) = o 3 = 1,...,N. (2.38)
Equation (2.37) says that linear momentum can only be lost due to interactions
with other phases. Equation (2.38) says that the interface can hold no linear
momentum.
2.3.3 General Case
Equation (2.1) is the constituent, microscopic general field equation. We
will find it useful to average this equation up to the mesoscale. The details
of the averaging are much the same as they are for the conservation of linear
momentum, so we will not repeat them here. Let us just say that the general
22
case, after averaging from the microscale to the mesoscale, can be written as
^(sQpQ^j) V (sQ{{ij)a + piva^ja pavH'ja))
+ V (Â£apaJV^a) ~ Â£apa~pa
)3
Â£aipJa f pp (waPj vj) na da
p? (wai)i
v^) na da
A
+ Â£apa>&a Â£apa^J ,
(2.39)
where the means that we are subtracting and adding this quantity as it will be
needed for certain variable definitions. We will find it convenient to substitute
the quantities from Table 2.1 directly into equation (2.39) as the averaging
procedure has already been performed.
2.3.4 Angular Momentum Balance
Angular momentum balance is probably the most difficult field equation
to upscale simply because the calculations are tedious. Thus, we will go into
more detail in this section than in others so that the reader may more easily
reproduce these results. We begin by substituting the appropriate terms for
the conservation of angular momentum from Table 2.1 into equation (2.39) to
obtain
23
Â£ da
+ eapa]r x (z + rV) Â£apajmJ .
(2.40)
Next we subtract x crossed with equation (2.27) (the conservation of linear
momentum at the microscale), where x is the macroscale field variable. To do
this we will need several identities which can be verified by combining (2.4) with
the fact that (2.6) (2.8) are integrals with respect to Â£ only, and not with
respect to x:
r x vi = x X vaj + Â£ x ib ,
vi{r x ib ) = x x vivia + tp(Â£ x tb) and
r x gJ = x x gai + Â£ x
among others similar to these. Starting with the left hand side, for the terms
we have
d
d
d
XV3) + (Â£>*Â£ Xvj ) XX (Â£apaJVai) = {Â£apajÂ£ XVl ).
(2.41)
For the V terms we have:
V [spa'vajx x v3 + Â£apaji>QjÂ£ x Â£ax x (tj)a
 Â£Q{Â£ X tj)a Â£apa]Va'X X V3 Â£/9jUQjÂ£ X V
+ sapajx x viyi + Â£ap3v^ x vi]
+ x X V [eq (Â£*') + Â£apa>VaiVa3 Â£ap3vh^ Â£apaiva3va3]
= V (Â£apa3Va3X X V*3) X X V (Â£apa3Va3Vaj)
24
 V (sax X (tJ)a) + X x V (Â£a{tJ)a) V (Â£q(Â£ X tJ)a)
V (sapaiva3x X vaj) + X X V (sapa,va,vaj)
4 V (eapajx x vjvia) x x V (sapajvivja)
+ V (Eapa*viÂ£ x via). (2.42)
The first and third lines after the equality sign in equation (2.42) cancel. Note
that if we switch to indicial notation
V (sax x (t3)Q) + a; x V (ea(tj)a)
(^ %itjkÂ£ijl\k T ijk\k^ijl
Â£ %i,ktjkÂ£ijl ^ ^i^jk^k^ijl T & ^i^jk^k^ijl
Â£ ^ik^jk^ijl ^ (2.43)
where the last equivalence is left in indicial notation because there is no equiv
alent in direct notation. We also have
V (Â£apajx x vJvJ ) x X V (eapQj ibib )
= (Sap:iXiVjVkÂ£lJi)ik Xi(sapajVjVk)jÂ£ijl
= Â£apa3Xi%kVjVkeijl + Â£ap3xi{vjvk)tkÂ£iji
 Eap3Xi(vJVk)jÂ£iji
= ^p^SikVjVkEijt = 0, (2.44)
and finally
V (Â£apa3vai(x x i>>)) x x V (Â£apaiva*vai)
= (Â£ap0,jVlXiVjÂ£ijk)tl Xi(Â£apaiViVj)jÂ£ijk
25
= Â£apa}XijViVjÂ£t]k + Â£apaiXi{ViVj)tlÂ£ijk Â£a p] X^Vj) j6ijk
P &ilVlVjÂ£ijk 0.
(2.45)
Using identities (2.43) (2.45) and treating the external supply of momentum
term, r x gi similarly, the left hand side of equations (2.40) becomes
d
dt
Â£apa*Â£ xvi ) + V {Â£apaiVj{$, x vi) Â£a(t X t3)a)
 Â£atjiSiji sapa>Â£ x gi .
(2.46)
Turning our attention to the right hand side of equation (2.40), and exam
ining the terms with no sum, we have
Â£apajx x (i + r3vJ ) Â£apajm3 + Â£QpQÂ£ x (i3 + r3vj)
x x Â£apaj (i r + vi )
= eapolm3 +Â£apa^x(i +?U) . (2.47)
Now, examining the terms on the right hand side which involve a sum, we have
\x x t3 + p3(x x v3)(wa0:i v3)] na da
[Â£ x t3 + p3^ x v3)(wa!3J ir7)] na da
[t3 + p?v3{wa0i v3) na da\
E
+ E
3^a
[\SVa\JA
Â£
xx J2
f3^c
Â£'
\SVQi JA
Â£
= E
0^a
[\sva
Wi Ja
 ( [f x tJ + pp(Â£ x v3)(wal3> ir7)] na da
(2.48)
26
Combining the left and right hand side, we now have the conservation of angular
momentum at the mesoscale:
d_
dt
(.Â£apa*Â£ Xl)J ) + V (Â£apQJVlÂ£ X Vl Sa(Z X t3))
Â£aI X (t3)a Â£apajZ X gj
= 8apajrnJ +Â£apax (i1 +?JvJ)
+ I 5./ I / [Â£ X t3 + ft (Â£ x v])(waf3i td)l nQ da
TZ LTO JA
8^c*
Equation (2.49) simplifies to
Â£a^Â£tjl = maj + M"' +
0jia
(2.49)
(2.50)
where we have replaced (t3)a by t> (we can do so because of equation (2.28))
and we have made the following definitions:
fnJ = Â£apajrnaj
(2.51)
is the rate of gain of angular momentum due to interaction with other species
within the same phase,
d
Maj = {eapai$, x vi ) V (Â£apajvj$, x v3 ) + V eQ(Â£ x t3)
+ Â£apaÂ£ x gi + Â£afPÂ£ x (i +Vvi)
(2.52)
is the rate of angular momentum gain by constituent j in phase a due to the
microscale angular momentum terms, and
= TJTT, [ [?xtJ+/f(^x v3)(wa0i v3)\ na da (2.53)
I^Kal JA
27
is the rate of angular momentum gain by constituent j in phase a due to inter
action with phase (5.
Summing over j 1,... TV yields the conservation of angular momentum
for the bulk phase
= Ma where the following restrictions apply: 3^ a (2.54)
N
tM 3) 9 II O Va, (2.55)
J2 m7 =0 a,/3^a J i v. (2.56)
Equation (2.54) implies that, in general, the mesoscale stress tensor ta is not
symmetric for multiphase systems. However, if iQj is symmetric then the right
hand side of (2.54) must be zero. Restriction (2.55) says that bulkphase angular
momentum can only be lost due to transfer to other phases, restriction (2.56)
holds because we assume that the interface is massless.
2.3.5 Conservation of Energy
Beginning as we did for the conservation of angular momentum, by substi
tuting the appropriate terms for the conservation of energy from Table 2.1 into
equation (2.39), we obtain
^(eapaj(ei + v3)) V (Â£a((t3 v3 + q3)a + pajva:i(eJa + vi)
1 _Q i__________Q
piv:>(ei + vi vJ) )) + V (eap,i;Qj(eJ + vJ v^))
Â£apaj(gj vi + IP )
t3 v3 + q3 f p*(e3 + yo3 v3 )(waf3] v3)
n da
28
1 /
=f(el HtP tr7) / pJ(wa'8^ rd) nQ da
2 J A
___ _______=a
+Â£apa,(Qj + V vi +?j(ej + vjvj) ). (257)
Now, from equation (2.57). subtracting equation (2.29) times v'1' and equation
(2.12) times (el + rPa), performing massive algebraic manipulations, and
again using equation (2.13), we obtain the conservation of energy equation at
the mesoscale:
Â£apai^(e*) V {eaqa*) Â£ata> :  Â£apaha>
= Qa+Qpj, (2.58)
where the colon indicates the tensor dot product (a : b = Yhij aijt>ij) Here we
have made the following definitions:
ha> = hJ + . v<*i (2.59)
is the external supply of energy.
a 1
I
2
e7 = e3 HvJ iPvaj vaj
(2.60)
is the energy density,
qai = (qj)a + (tJ vj) ta va' + pa>va(ea> + va^ )
pajv3(ei + vi vi)
(2.61)
is the partial heat flux vector for the jth component of phase a
^  Q ' l ^
Qaj = gapCij^Qj __ p VJ ( ^ _ipvj vQi^P ) VaJ Iri(el + vi vi)
1
r3 (eJ + vi vJ ))
(2.62)
29
is the rate of energy gain due to interaction with other species within the same
phase not due to mass or momentum transfer, and
J A
vi)(walii vJ) na da
vaj / [i7 + p?vJ(wal3J v7)] na da
(2.63)
A
is the rate of energy transfer from phase (3 to phase a not due to mass or
momentum transfer.
The bulk phase conservation of energy equation is found in an analogous
manner to the bulk phase conservation of linear momentum by using (2.20) and
subtracting va times (2.33) and (e + v) times (2.21) and then summing
over constituents. This yields
Cfc fj f LX
The first restriction says that bulkphase energy can only be lost due to transfer
with other phases, and the second restriction states that the interface retains no
energy.
Pyt a
where the following restrictions hold:
= 0 j = l,...,JV.(2.65)
0 Vq, and (2.64)
30
2.3.6 Entropy Balance
After substituting the appropriate quantities from Table 2.1 into (2.39) the
entropy balance at the mesoscale for constituent j in phase a is
 v (sa((va*rpc
 paiviip )) + V (Â£pnJvnhf ) Â£apaib>
= SÂ£Pai  / & + fPrf{waPi vJ) na da
0^a ^ ' a' JA
T , f p>(wapi v>) na
P^WJa ________________
+ Â£a paj~KJa + rp + r'vj'.
da
(2.66)
Subtracting rf* times (2.12) and using (2.13) we obtain
eapa> V ' ~ Â£Paibaj
(267)
0
where we have made the following definitions:
Aa^=7Ja (2.68)
is the entropy production,
/v 
= Tp
is the entropy of the jth constituent in the aphase,
(j>a = (4*>)a + pv0'3rf] pajvjrp
(2.69)
(2.70)
31
is the partial entropy flux vector for the jth component of the aphase,
b>a = ba* (2.71)
is the external entropy source,
$3 = 77777 [ [& + ft ft J vJ)] n da
77777 f ft(wa!3i vj) na da (2.72)
0'/a Ja
is the entropy transfer to the jth component in the aphase through mechanical
interactions with the same component in the ,5phase, and
rf=Â£apaj(r}ia+ rr)j r3 rjaj) (2.73)
is the net entropy gain due to interaction with other species within the same
phase
The bulkphase counterpart is given by
Â£apan*r ~ v'{Â£ar) Â£apab = E ^ + AQ> (274)
where the following restrictions apply:
N
+901 V3) = 0 Va, and (2.75)
j=1
2 (8;j + e^VJ) =0 3 = 1 (2.76)
q,/3^q
The first restriction states that entropy can only be lost due to interactions with
another phase, and the second restriction states that the interface can hold no
entropy.
32
2.4 Standard Assumptions and Resulting Entropy Inequality
In this section we make a few standard simplifying assumptions, and give the
resulting form of the entropy inequality. Henceforth we assume that the stress
tensor, iQ, is symmetric. We also assume that the system is in local thermal
equilibrium. That is, because the phases are viewed as overlaying continua, we
assume the temperature of all constituents in all phases at a single point is the
same for all phases. This assumption can be expressed as
T = Ta*=Tfii Vq, V/3,Vj. (2.77)
Note that this does not mean that we are assuming that the temperature is
constant; T is still a function of time and space.
Next, we assume that the system is thermodynamically simple in the sense
of Eringen [24], This means that entropy flux and external supplies of entropy
are due to heat flux and external supplies of heat alone, respectively. This
assumption can be expressed as:
= (2.78)
hai
bai = . (2.79)
Entropy is a mathematically useful quantity. However, experimentally it cannot
be measured directly. Thus, we choose to perform a Legendre transformation to
convert the internal energy ea* into the Helmholtz free energy An,
Aa =ea> +Trf, (2.80)
33
allowing us to choose temperature instead of entropy as an independent variable.
The second law of thermodynamics states that the total entropy generated
by the system must be nonnegative, and is maximum only when the system is
in equilibrium. This statement can be expressed as:
a = A" = Aj 
a j1
(2.81)
We begin by solving equation (2.67) for subtract A times the con
servation of energy equation (2.58), to eliminate the heat source variables, hai,
perform the Legendre transformation with (2.80), then sum over all constituents
and phases to obtain
sapa (DaAa aDaT
Dt
+ r
Dt
V j=l
1
"f
r N tai vaJQ 
l 3=1
} + T + V(Â£nf/^An')
3=1
+ Y (*J pajAajI) :
3 1
1 v> tv
j=i
B*a V ' )
(2.82)
34
where vaja = va, and in general a comma in the superscript denotes a
difference in the superscripted quantity. Identities needed to obtain equation
(2.82) can be found in Appendix A.3.
2.5 Choosing Constitutive Independent Variables and Exploiting
the Entropy Inequality
Thus far, the theory introduced applies to all media meeting the require
ments laid out in the beginning of this chapter and the previous section. What
distinguishes one material from another is the set of variables upon which all
other variables depend; these are called constitutive independent variables, or
more simply, independent variables. How we choose the independent variables
is informed by knowledge of the material we wish to model, experience, and
experiments. It is important to note that the choice of independent variables
defines the material.
Once the constitutive independent variables are chosen, we expand the
Helmholtz potential term found in (2.82) in terms of these variables, perform
various algebraic manipulations and collect terms. The axiom of equipresence
requires that initially all constitutive variables be functions of the complete set
of independent variables, even, for example, if one phase lacks a certain con
stituent. Only after the entropy inequality is exploited can the concentration of
a species be set to zero, [31]. Failure to adhere to the axiom of equipresence can
produce erroneous results. However, it can be shown that the Helmholtz free
energy is not a function of all the independent variables, as will be demonstrated
momentarily.
35
To illustrate how we chose the independent variables and exploit the entropy
inequality in the sense of Coleman and Noll, [19], we will present two illustrative
examples. The first example will be for a singlephase elastic solid composed of
a single constituent. The second example will be a singlephase thermoviscous
fluid, also composed of a single constituent.
2.5.1 Example: Elastic Solid
Mechanics of materials tells us that linear elastic bodies exhibit different
behavior depending on the strain and temperature of the material, denoted E8
and T, respectively. Expanding the derivative of the Helmholtz potential as a
function of these variables we have:
DAS dA8 DES dAs DT
~Dt ~ dE^ldf + dTDF ^
where we have dropped the superscript on the derivative terms since there is
only one phase. Substituting this into entropy inequality (2.82) and collecting
like terms we obtain:
1 ( cb4s ,s /j?s\T\ gs
T dE
 {F8)1 ts (Fsy
VT
+ T2^S> 0.
(2.84)
where we have converted ds into the independent variable Es via the identity
ds = (Fs)~t if (Fs)1.
Variables which appear linearly in the entropy inequality which are neither
constitutive nor independent can vary independently, thus their coefficients must
36
be set to zero to avoid violating the entropy inequality. In this example, the
following variables satisfy the these requirements:
T, E\ VT.
(2.85)
The resulting equations are called nonequilibnum results as they must hold
both at equilibrium as well as away from equilibrium. Since the coefficient of T
must be zero, we obtain the classical result
dAs
Iff
I ,
(2.86)
which means that temperature and entropy are dual variables. For the coefficient
of if, we obtain
dAs
t$ = pSFS dE* ^ (2'87)
which is a classical Hookes law. Finally, for the coefficient of VT we obtain the
result that
qs = 0. (2.88)
This last result says that there is no heat flux. As such, the material that we
have defined by our choice of independent variables does not conduct heat. This
is perfectly acceptable if we mean to model an insulating solid or we know a
priori that the system has a constant temperature. However, we may wish to
model an elastic solid that conducts heat, that is, a thermoelastic solid. To
accomplish this, we add VT to our list of independent variables, resulting in an
37
additional term in entropy inequality. Equation (2.84) becomes
 (FT1 f (FTr) : E'
, ps dAs D(VT)
+g ~Td(VT) Dt ' ^2'89)
The variable VT no longer appears in list (2.85) because now it is independent.
In its place D^tT^ appears, leading to the nonequilibrium result that
dAs
Ps ( ' dAs
T l ^ dT
1 ( J)Al
t{ 1 dE
VT
qs 
= o,
(2.90)
0(VT)
meaning that the Helmholtz free energy is not a function of VT. Where it
can be demonstrated that Helmholtz free energy is not a function of certain
variables, as above, the Helmholtz free energy will be taken as a function of
a subset of the complete list of independent variables to simplify calculations,
otherwise we will adhere to the axiom of equipresence.
We now have a heat conducting elastic solid, that is, a thermoelastic solid.
Equilibrium for this systems is guaranteed to occur when VT = 0. At equi
librium the net generation of entropy is minimum, so that 0 and
a(Vr)2 Taking the partial derivative of equation (2.89) with respect to
VT and setting it equal to zero yields
9s = 0,
(2.91)
that is, there is no heat flux at equilibrium, which makes sense for a thermo
elastic body.
The coefficients of the variables that are zero at equilibrium are a function of
these variables. To clarify, if x and y are variables becoming zero at equilibrium,
38
then the coefficients of x and y are functions of x and y. One way of obtaining
nearequilibrium results is to linearize the coefficients to form positive quadratic
terms by using a Taylor series expansion about the variables becoming zero at
equilibrium, and then truncating all secondorder and higher terms. In this
simple example, we have that qs = qs(VT), so that
q = q%q + Ks VT, (2,92)
where eq stands for equilibrium, and Ks is a secondorder tensor resulting from
the linearization process. Equation (2.91) tells us that qs\eq = 0, thus
qs = Ks VT. (2.93)
Equation (2.93) is Fouriers law of heat conduction, which is known to hold for
thermoelastic solids near equilibrium.
2.5.2 Example: Thermoviscoelastic Fluid
If we want to model a thermoviscous fluid we would include in the list of
independent variables the density, the temperature, the gradient of the velocity,
and the gradient of the temperature, denoted pl, T, Vvl, and VT. respectively.
However, it can be shown that Vvl is not frame invariant, where as the sym
metric part of Vvl, d! = l/2(Vvl + (vl)T) is frame invariant, [24]. Thus, we
include d! in our list of independent variables instead. Including gradient of
various variables produces flow in those variables. For example, just as includ
ing VT produces heat flow, including a gradient in the volume fraction, Vs1,
produces a material which has flow dependent on the moisture content, which
is the case in swelling materials. Including higherorder gradients of a variable
39
Using the method outlined in the previous example, we obtain the following
nonequilibrium results
dAl
~dT
(2.98)
which states that entropy and temperature are dual variables, and
dA1
ddl
= o,
(2.99)
which means that the Helmholtz free energy is not a function of d!.
We define equilibrium to be when both dl and VT are zero. If a: is a variable
becoming zero at equilibrium, then = 0 because the net generation of entropy
of the system must be minimum at equilibrium. Thus, we obtain the following
equilibrium results. For the coefficient of dl we have:
tl = p'l, (2.100)
where we have used the thermodynamic definition of pressure
p = (2.ioi)
Pressures will be discussed at great length in subsequent chapters. Equation
(2.100) is a standard result and states that at equilibrium the stress in the
liquid is equal to minus the pressure in the liquid phase. For the coefficient of
VT we have
ql = 0, (2.102)
which states that at equilibrium the heat flux is zero. Since we now have
more than one variable defining equilibrium, the condition on the second par
tial derivative of A in the previous example is replaced by the requirement that
41
the Jacobian of A be positive definite. This assures that entropy production is
minimum.
Once again, we would like to obtain nearequilibrium results for the system.
The situation is different than the previous example because now we have more
than one variable that defines equilibrium. We will still use a Taylor series
expansion about the variables which define equilibrium and truncate all second
order and higher terms. However, there is now more than one way to accomplish
this. First, we can choose to do single variable expansions. Setting the coefficient
of d! to Ql, we have that ql ql( VT) and Ql = ql{dl), analogous to the previous
example. Thus, single variable expansions yield
ql = KlVT (2.103)
tl = plI + Ll : d!, (2.104)
where ql\eq = 0 because of equation (2.102), and Ql\eq = 0 because of equation
(2.100). In the above equations Kl and i) are secondorder and fourthorder
tensors, respectively, resulting from the linearization process. We note that
because we choose to do single variable expansions that Kl and Ll must be
evaluated as functions of dl and VT, respectively. In general, whenever single
variable expansions are performed the coefficients resulting from the lineariza
tion process must be evaluated as functions of all the other variables which
define equilibrium. Of course, having constant coefficients is preferable and can
be achieved by using expansions in all of the variables that define equilibrium. In
this example we can perform two variable expansions by taking ql = ql(VT, d!)
42
and Ql Ql{dl, VT). Doing so yields
(2.105)
(2.106)
ql = Kl VT + H1 :dl.
tl = ftI + Ll :dl + Ml VT,
where Kl and Ml are secondorder tensors, and Hl and Ml are fourthorder
tensors resulting from equilibrium expansion, all of which are only functions of
independent variables not becoming zero at equilibrium. Unfortunately, using
Taylor series expansions in all the variable defining equilibrium can produce
cumbersome expressions. We will be careful to the use expansions that are
relevant to the system we want to model but significantly enrich the resulting
equations.
In subsequent chapters we will use the techniques described in this section to
obtain nonequilibrium, equilibrium, and nearequilibrium results for far more
complex systems.
43
3. A New Choice of Independent Variables
Deriving a physically meaningful transport equation for porous swelling ma
terials that undergo finite deformations depends largely on our ability to relate
thermodynamically defined variables to physically interpretable quantities. In
this chapter we present a novel and judicious choice of independent variables for
the solid phase that clarifies the relationship between thermodynamically de
fined pressure and actual physical stress. To further show the usefulness of this
new set of independent variables, we will show how it elucidates and simplifies
the assumptions and the derivation of the transport equation previously inves
tigated by Singh et al., [51], while continuing to capture the important features.
We then show how this framework can be used to derive a transport equation
used to model swelling polymers, and compare the result with another transport
equation derived using the FloryHuggins theory.
3.1 Motivation
In this chapter we will restrict our discussion to a twophase system consist
ing of a viscoelastic solid and viscous liquid phase, each composed of a single
constituent, and neglect all interfacial effects, though the theory can be used to
incorporate such effects, [7, 8]. Results for nonswelling media may be obtained
by making simplifying assumptions.
To capture the viscoelastic nature of the solid phase and phase interaction
effects, the independent variables should include variables upon which the solid
phase behavior depends, the liquid phase behavior depends, and the interaction
44
of phases depend. The intuitive choice for independent solid variables include
but are not limited to: the solid phase volume fraction, es, solid phase density,
ps, and smearedout solid strain, Es. These variables are not independent. To
see this, consider the continuity equation for the solid phase assuming no loss
of solid material to the liquid phase, see equation (2.23):
Ds(ssps)
^+Â£yV.v=0. (3.1)
Let Fs denote the deformation gradient. Then the right CauchyGreen tensor
can be expressed as
Cs = (Fs)t Fs, (3.2)
where Fs is the deformation gradient of the solid phase and is defined by (in
indicial notation)
Fs = FskK = x
Lk,K
dxl
ox
(3.3)
K
where xs is the Eulerian coordinate and Xs is the Lagrangian coordinate, and
where the strain tensor is expressed as
Es = 1/2{Cs I).
(3.4)
Let V be the total volume of the smearedout solid phase, and let Vo be the
initial volume of the (smearedout) solid phase. Equation (3.1) can be expressed
in integral form as [14]
/,
espsdv
e3psJsdvo
(3.5)
45
where the jacobian, Js = det (FsT Fs) = det(Cs) = det(2E3 + I). From
(3.5), we have J3 = Combining these results, we obtain the following
relationship
det(2Â£s + I) = Js = (3.6)
Â£SpS
Thus, we see clearly that E3, p3 and Â£s are dependent. In particular, the first
equality tells us that all six components of 2Es + I are not independent of Js.
Because of this coupling it is not straight forward as to how to choose inde
pendent variables which lead to physically meaningful variables. A first choice
might be to eliminate volume fraction or density as one of the independent vari
ables. However, density is dual to the thermodynamic definition of pressure, and
volume fraction measures liquid content, and hence both are essential quantities
in determining macroscopic behavior.
One way of handling this coupling is to enforce the continuity equation using
a Lagrange multiplier [33]. The Lagrange multiplier then becomes an unknown
of the system to be determined by boundary conditions [24]. Using a Lagrange
multiplier, derivatives with respect to Es (or Cs) can be evaluated letting p3 and
es vary without restriction. However, it then becomes unclear what the physical
representation of the dual variables are. To see this, consider the following result
obtained via HMT using ps, Â£s, and Es in the list of independent variables [49]:
Fl
= psI + tse + tsh, (3.7)
where
r)As
t3e = p3F3~(F3)T, (3.8)
46
and
(3.9)
In equation (3.7), ts is the stress of the solid phase and is related to the stress
and tse and tsh are labeled the Terzaghi stress and the hydration stress, respec
tively [11, 12]. As 1/3tr(ts) represents the physical pressure exerted on the
solid phase it is not clear what portion of pressure ps represents, since neither
tse nor tsh have zero trace.
Another way of dealing with this coupling is to make simplifying assump
tions. For example, in [49], Singh et al. proposed a constitutive theory to model
a twophase polymeric system with a viscoelastic solid phase and viscous fluid
phase. By including strain, time derivative of strain, and a gradient in strain as
constitutive independent variables, they developed a theory applicable to vis
coelastic systems. Therein, they obtained novel forms of Darcys law. In [50],
they extended the twoscale theory to three scales, and then in [51], a trans
port theory was developed based on their work in [50]. The development of the
transport model required the assumption that the strain of the solid phase and
rates of strain, were a function of volumetric changes only, that is, a function of
Js, and not a function of shear.
This brings us to another issue we address in this chapter. The form of
Darcys law used to derive the transport equation in [51], which we will derive
of the entire medium via t = ests + eltl. The variable ps is the thermodynamic
pressure and is defined by
(3.10)
47
in a subsequent section, is written in terms of gradients in the classical pressure,
volume fraction, and time derivatives of volume fraction. However, as we will see
later, this form of Darcys law can be written in an equivalent form expressed in
terms of a gradient in the Gibbs free energy (or chemical potential) of the liquid
phase. Either form of Darcy's law, when combined with the bulk conservation
of mass equation for the liquid and solid phases, yields a transport equation.
However, we need different constitutive equations and employ different assump
tions depending on the form of Darcys law that we use. We will derive two
such transport equations and discuss the similarities, differences, and usefulness
of these equations.
In this chapter we will first introduce a new set of constitutive independent
variables that give rise to a physically meaningful interpretation of the solid
phase stress, solid pressure, and Terzaghi stress tensors. We will present the
necessary techniques and methodology to obtain the entropy inequality and then
present only the novel results stemming from this choice of variables. Then,
based on this choice of independent variables we derive a transport equation
analogous to that found in [51]. Next, using an equivalent form of Darcys
law, we will use FloryHuggins theory to develop another transport equation
specifically applicable to polymers. Lastly, it is worth mentioning that we are
able to reproduce the three pressure relationship of Bennethum et al. [14] in
terms of these new variables.
3.2 Constitution
As discussed above, complications arise from the dependence of the vari
ables Â£s, ps, and Es through the solid phase continuity equation. In choosing
48
an appropriate set of independent variables consider the properties of swelling
clays and polymeric materials. Swelling clays, as well as many polymeric mate
rials, have markedly different behavior at low and high moisture contents. At
high moisture contents the liquid phase cannot support shear. Thus shearing
has little affect on flow, and macroscopic behavior is dominated by volumetric
changes due to changes in moisture content. In contrast, a swelling clay with
low moisture content will have very few layers of vicinal water (water residing
close to the solid phase), and the liquid phase supports shear and the effects
of shearing become of greater concern when modeling deformation and trans
port [23]. For these reasons, it becomes appropriate to consider splitting the
deformation of the (smearedout) solid phase into volumetric (dilatational) and
isochoric (distortional) parts. Ideally it would be nice to decompose strain, Es,
into volumetric and isochoric components, however mathematically this makes
the problem quite complicated due to the relationship between es, ps, and Es.
Consequently, we follow Holzaphfel [32], and consider a multiplicative decom
position of F3 and Cs as follows:
Cs = (JS)2/3CS, (3.11)
Fs = (Js)1/3Fs , (3.12)
where Jsl^I and Js2//3J represent volumetric deformation, Fs and C$ are the
modified deformation gradient and the modified right CauchyGreen tensor, re
spectively. These variables account for distortional deformation. They are re
lated in the same way Fs and Cs are, that is,
Cs = (Fs)tF3. (3.13)
49
With this decomposition in mind, we assume the following independent vari
ables:
sl, Vs1, (e', V(sl, T, VT,
___ Jn)
pl, Cl>, cs, cs Js, Cs\ Va's, via, df, (3.14)
where m = 1,... ,p and n = 1,... ,q denote material time derivatives of order
to and order n, and j = 1,..., N 1 represent N 1 constituents throughout
this thesis unless otherwise noted. This choice of independent variables differs
from that found in [49] in that we have replaced Es and ps with Cs and Js.
Additionally, we have replaced higher order derivatives of Es with derivatives
of el and Cs.
3.3 The Entropy Inequality
In this section we will present the necessary methodology to obtain a form
of the entropy inequality that is exploitable.
To simplify the following manipulations, we postulate the dependence of the
Helmholtz free energy as follows:
A1 = A1 (sl,
. , _ M
A* = As(sl, il, JS,CS,T, C\CS),
(3.15)
(3.16)
where to = 1,... ,p and n = 1,... ,q are material time derivatives of order p
and q, Caj is the concentration of the jth species in phase a and j = 1,..., N.
Otherwise, we adhere to the axiom of equipresence, that is, all other variables
are considered a function of the complete list (3.14). Material time derivatives
50
of A1 and As appear in the entropy inequality (2.82). Using the chain rule, they
may be calculated as follows:
DlA dAl D'e
y ^ dAi
dAl Dl 'e
3Al Dl
+
Dt ds1 Dt Dt dp1 Dt
m= 1 Uc '
Nl Dl(Cl>) Dt
dAlDlT 8Al DlCs A c _ _ .
+  + ^ + Â£ gr : Di.
n=i q Cs
<)T Dt dCs
Dt
(3.17)
DM* dAs Dsel dAs Db{s)i dAs DSJS DS(CSD
+ Ers;Dr + 8jrDr + E''J
m= 1 U C
Dt
ds1 Dt
dAs DST dAs DSC
+
&T Dt + dC
3 = 1
Jjl)
q dAs Ds Cs
Dt
Dt
Eosi
^77
Dt
(3.18)
"=1 d C
Because j = 1,...,JV in (3.15) and (3.16), the concentrations, CJ, are a depen
dent variables since
N
= (3.19)
3 = 1
which implies for example
DaCaN
Dt
Nl
E
3=1
DaCaJ
Dt
(3.20)
giving rise to the relative chemical potential [13]
= pa> /zQ*,
(3.21)
where
A 3 =
dAa
dCai
(3.22)
We will find it useful to be able to convert a material time derivative with
respect to one phase (or species within a phase) into a material time derivative
51
with respect to another phase (or species within a phase). This can be done
using either of the following two identities:
Da() Ds()
or ~isr+*'*<) <323)
Â£>">() Da()
m=m+v"v{) <3'24>
Since ds is not in the list of independent variables (3.14), it is necessary
S
to convert dx into the independent variables C and Js. We use the following
calculation:
ds = (ftt
i
Es (F'
S\ 1
= 1(FTTCS(F)1
= \(FS)~T (^{JS)1/3JSCS + {Js)2/3&^j (F8)"1
= J(Js)1>(Fs)t cs (F8)1 + \(FS)~T ts (F3)1
= + (fs)~t c3 (Fsy\
(3.25)
In the above calculation the first equality is an identity, see [37], the second
equality is due to the definition of the strain tensor, the third and fourth equal
ities are due to equations (3.11) and (3.12), respectively, and the fifth equality
holds because of (3.2).
Hence it follows from (3.25) that
= \ ((Fr (ftT) # + Â£ <326)
52
and
: dS = \ (^r1 EE : s
N
+ ;^Ewfo.
(3.27)
j=i
Using the identities given in the entropy inequality (2.82) above and collecting
terms, we obtain the following form of the entropy inequality:
DST
1
"f
f
/ch4a \ U +E
s
N Ic< T"
lE^ /A'C'1 
Dt
a j=1
_I lrsjdAa' ^
T\PdJs 3
1
'T
i id A1
Â£pl =s + Â£ p
dC dC
T dp1
Etr(i^)) j
j=i /
8AS Â£s ,s
N
(Fy
E tsn(Fsy
\j=i /
s
: C
A { *
+f (Â£1:d'
\j = 1
tÂ£(x>v^Â§)
a \n=l Q(f)
a \ (+i)
:CS
1
T
, (dAl A , tdA1 t dAl
Â£ p \W + rilVT + Â£ p 7^VÂ£ + Â£ : V C
del
n=0
dc
+*V Â£ h v1?' + t'p'gTW+Â£ VVC + t[
m=i del U/J j=i
Z,3
53
1 f jjdA1 ^s^sdAs \ (m+D,
m=1
rE('V^+^V^
N
\a Â£aVT I v>
+X!^~ {q + J2
l j=i
paJVa*'a ( AQ + {va>'af ) taJ Va>'a
fEE [r?+r + v(Ev^
Q j=l
+EE^ ^QjJ): (v^q)
Q j=l
+f+
V
J=1
f Eis(4n
0^a
(3.28)
We choose to use Lagrange multipliers to weakly enforce the continuity
equations, see [33]. The entropy inequality is modified as follows:
AQ
Anew A0ld + ^ rj,
Da{Â£apa)
Dt
+ Â£apaVVa
a,a^/3
N
+
EE
a j=1
\ai
~Â¥
DCai
aa + V (eapa*va>'a)
Â£ P
Y eg + r3 ~ Ca> Y ?
\a,Ck^0 Ct,Ck^fi
'a
/3
> o.
(3.29)
where Aqm is given by (3.28). This yields the following form of the entropy
inequality:
A =
_1
T
a
8Al
dA
dT
t t^rr s S9A3
Â£P^~Â£P^
+ 77'
D'T
Dt
 pXl  r
54
~ E ^E &*'
1 ( I idA1 A
t(Â£pW~Â£X)p
1
T
1
T
. <)AS 1 Â£
N
Â£spswlr
dJs 3 Js v >
tj^4'
j=i
. s^s
+ Â£sps=s
Js
ac &c
l
~i ( n \
r ] : d1
Â£sps Xs
(Fsrl (Fsy
a
: C
+ T
^.7=1
tE(x>v"^)
V"1 ac'/
dA \ ^
I u
1
T
e'',(^ + ')VI'+Â£V\
^ ^ V,
JL f)Al (") P_ , .
/ / X ^ t'ZT __, p=rS / / X "V t/yi __(m)i
+Â£V Ey:VC +Â£V E^Ve
n=() ^ Â£s m 1 C' Â£
+
I V (Vp1 +
TÂ£iv%y yv
vls
v^eQVT v
+Ew E
L j=i
9 Â£ 1
iV r
rp2
a
,a Af1
+ EtE
Or J = 1
Â£> vQj,a pajvaj'a ( + (i)aj,)2
p'
r>   p"'(/r> yiQjv ) + pn(\a' aow)
: Va'"
?EEi<*? +"J) ^<7 + ?*>+w
Nl
o j = l
f)
\paNJ
55
(AQ^ AQiV)V(fQpaO + V [EapQi{Aa> Aar
1 N
4d
J=1
Aft' +
DELES'
Aaj + (uJ,Q)2
2V 7
iV1
q j=l (3^a
~^EE^ ,4a + ^K'!)2 + AgrC'
q L j=1
where we have used the following two identities
N N1 ,
y^ paj. vaia = y^ / paj Â£_Lp<*N j.
t=i j=i ' p A
> 0, (3.30)
(3.31)
TV1
1=1
3~ 1
GQ^ : VvQ^a = E ( G' ~ ^GN) : Vv
N1 / a \
Gaw Vv i7aiQ,
VPQAV
(3.32)
to remove the 1V
[13].
3.3.1 Two Phase Pressures
In Gibbsian thermodynamics, pressure is defined as follows:
8 An
r
dV
(3.33)
Mai...
where Ma> is the mass of constituent j in the aphase, Va is the volume of the
aphase, and A? is the total (extensive) Helmholtz free energy. Dividing the
top and bottom of the partial by the magnitude of the REV, we obtain
dAa
_a __ _d{eapaAa)
dea
= Â£apa
de
Â£p
(3.34)
ep
56
so that we now have an expression for the thermodynamic pressure in terms of
intensive variables. Another thermodynamic variable with arises is
j=i
Bennethum and Weinstein showed that [14]
pa = pa
(3.35)
(3.36)
for a singlephase material made up of one constituent by converting to extensive
variables. However, as we will see, these pressures differ for a swelling porous
media. Exploitation of the entropy inequality in the sense of Coleman and Noll,
[19] yields the following equilibrium result for the coefficient of dl:
tl = plI, (3.37)
which implies that
Itr(t') = pl. (3.38)
The stress tensor is a physical quantity, and since one third the trace of the stress
tensor is the actual physical pressure exerted on an isotropic system, equation
(3.38) indicates that pl coincides with the physical pressure in the liquid phase.
Thus, it is called the classical pressure (see also [14]).
The relationship between the thermodynamic pressure and classical pressure
was elucidated in [14] by considering the Helmholtz potential to be a function
of two different combinations of the same set of independent variables. As such,
they are not different potentials; an overline is used to distinguish between the
57
sets of variables being used and will be dropped later. Consider the following
relationship for the liquid phase:
A'ie1 ^ ,...) = Al(el,pl,Cl>,...),
(3.39)
where j = 1,..., N for Â£lplJ and j 1,..., ./V 1 for Cli, with the total differ
ential given by
dAl
DA=i
N
, ^ 0Al
, h Â£+p;d(zlplj)
Elp
d(Â£lph)
dA
del
p,ci
, dA
dÂ£+w
, i ^ dA1
^ +E och
r'.C'j,... 3 = 1
dC1* (3.40)
Then, for example
rl
dA
dsl
dA1
del
dA1
del
N
dA1
. d(Â£lfJi)
j=i v '
dA'
sl d(Â£lph)
epli,... 3 = 1 V r 7
d(Â£lpli)
dÂ£l
(3.41)
Multiplying (3.41) through by Â£lpl we obtain the three pressure relationship in
[14]
pl(Â£l,pl,Cl) = pt(Â£l,Â£lpl>) + nl(Â£l,pl,C1),
(3.42)
where
t jid A1
7T = Â£ p
del
(3.43)
is called the swelling potential since it represents the degree to which the energy
of the phase is changed with respect to a change in the liquid content. Equation
(3.42) states that the classical pressure is equal to the pressure obtained by
58
changing the volume of the liquid phase keeping the mass fixed plus the pressure
obtained by changing the volume fraction keeping the density fixed.
In a similar fashion, one can show that
N
pl{Â£l,Â£lph) =
3 = 1
, ,.dA'
pp3w>
l\2
dA
dp1
= p(Â£l,pl,C,i), (3.44)
so that the definition of the classical pressure coincides regardless of which set
of independent variables is used.
Using similar calculations one can show that equations (3.42) and (3.43)
hold with l replaced by s, so that we obtain the general relationships
and
7T
apa
dAa
de'
a
l,s,
(3.45)
pa(ea,pa,Cai) = pa(Â£a,Â£apai) + na(Â£a,pa,Ca), a = l,s. (3.46)
The liquid phase pressures remain the same as in previous works [14], since
the independent variables have not been changed. However, in the solid phases
we have replaced the variable ps with Js. Consider, the following relationship
As(Â£l,ps,Cs>,...) = As(Â£l,Â£sps\...) = As(Â£l,Js,Csi,...), (3.47)
where j 1,..., IV 1 for CSj, and j = 1,..., N for Â£spS]. As before, the
overline and tilde are used to distinguish between the different combinations of
59
the same set independent variables. The total differential is given by
N1
DAS =
dAs
8sl
dAs
ds1 +
dAs
p.CaJ
8ps
dP
8sl
8Aa
N 8AS
ds1 + Y .
^ d(espsi)
epsj j=1 V r
,CJ J=1
d(sY)
dCs
8sl
, dA
ds ~t~
J,csi
8JS
N1
" + E
dAs
eCsJ
3 = 1
dCs>. (3.48)
Then, for example
8AS
8es
8 A
ps,CsJ
8es
N
pS] 8 A
Â£ap
^ ss 8ps
'j 3 = 1 F
8AS
8es
zl ,js
Js 8AS
ps,Csi
ss 8JS
,(3.49)
Â£\C33
where we have used the following calculation
8Ja Â£oPod{l/es)
ds3 8ss ps 8ss
ps,Ca3 p\CS3
Â£oPo Ja
pS,CsJ
(ss)2ps Ss
Multiplying equation (3.49) through by esps we obtain
,8AS
(3.50)
7T* = pS pS = pS + pS JS
8JS
(3.51)
,C*i
where the first equality is the three pressure relationship (3.46) with a = s, and
the second equality yields the following identity for the classical pressure of the
solid phase in terms of the new set of independent variables
,8AS
pa{sa,Ja,Ca>) = psJs
8JS
(3.52)
Â£\CS3
Using similar calculations one can show that
na{Ja,pa,Csi) = paJ!
8AS
8JS
(3.53)
p\C 3
60
where we note the difference in the variables being held fixed, and that
pS{Â£S,JS,CS>) = SSpS
c)As
des
(3.54)
3.4 Novel Results
Rather than present the entire constitutive theory associated with the ex
ploitation of entropy inequality (3.30) we will present only those results which
are novel and relevant to the current discussion. A comprehensive discussion
of these results can be found in [49]. We begin by addressing how the choice
of independent variables, list (3.14), affects the form of the solid phase stress.
Next we present two equivalent forms of Darcy's law. One form of Darcys law
will be used to derive a transport equation that is analogous to the transport
equation derived by Singh et al., [51], but where the new choice of independent
variables simplifies and elucidates the derivation and assumptions necessary to
obtain said equation. Using a second form of Darcys law, we will derive another
transport equation using FloryHuggins theory for polymer solutions. We end
this section by comparing and contrasting the resulting transport equations.
3.4.1 Solid Phase Stress
After exploiting the entropy inequality in the sense of Coleman and Noll,
this choice of independent variables produces the novel nonequilibrium result
for the coefficient of Js:
Itr(t') = p\ (3.55)
where we have assumed that the diffusive velocities in the solid phase are negli
gible and ps is given by (3.52). This result indicates that ps now coincides with
61
the classical pressure in the solid phase, whereas before it was unclear what
portion of the pressure ps represented (see equations (3.7)(3.9)).
a
We obtain the following equilibrium result for the coefficient of C :
^ dAl
= ps\sI + 2plFs ^ (Fsf + 2pF ^ (Fs)t, (3.56)
where Xs is a Lagrange multiplier used to enforce the continuity equation. Taking
onethird the trace of equation (3.56), using (3.55) to eliminate ts, and solving
for ps Xs we obtain
dAs
2 slpl dAl s 2 dAs
F ^ 3 Â£s 8C 3 dC
Substituting this equation back into (3.56) we have
Â£l rsh yse 2 Â£lpl 8A s 2 3 8AS
ts = psI + t*'1 + t
Â£'
3 es dCs
C 77 P:
3r dC
(3.57)
(3.58)
where
tse = 2 psFs d A3 dCs (^)r, (3.59)
tsh = 2plFs d A1 ocs (fs)t. (3.60)
While it is possible to lump the last two terms of (3.58) into the definitions of
tsh and tse so that they both have zero trace, we choose to keep with convention
and define equations (3.59) and (3.60) in an analogous manner as equations (3.8)
and (3.9). Note that now, when we take one third the trace of equation (3.58)
we recover equation (3.55).
62
3.4.2 Darcys Law
Darcys Law is obtained by combining the linear momentum equation for the
liquid phase, equation (2.36), with constitutive equations for the stress tensor
and the exchange of momentum term, Ts. A complete mesoscale nomencla
ture is given in Appendix A. Exploiting the entropy inequality results in the
equilibrium result
i
m= 1 O Â£ 1 n=0 Q ~Â£
where pl is defined by equation (3.34) with a l. Thus, linearizing about the
variable vl,s vl vs, we obtain
T = Rl (ielvl's) + pfVe1 s'p1 V V {Â£)l
3 \ / r r _(m).
in 1 0 Â£
JL <2+l>
 E % :
n=0 q Cs
To eliminate Ts and tl, equation (3.62) along with the nearequilibrium result
for tl
E
dAl
("+!)
: (V c y
(3.61)
I I I i ^> DA __(m)/
Ts = p'VÂ£lÂ£lplTVÂ£l
elP
tl = plI + L1 :dl + H1 VT
(3.63)
(see [49]), are substituted into the conservation of linear momentum equation
(3.61). The result is Darcys Law:
R1 (sV3) = V(eV) + plVsl + V (elLl : dl) + V (elHl VT)
9 dAl ,,^dAl.
i J.
+Â£ pg
<"*>/
n=0 Q Â£ts m= 1 D Â£
(3.64)
m= 1
63
where Rl is called the resistivity tensor and is second order. On the right hand
side of the equation the first term is responsible for flow due to a gradient in
pressure and is the primary driving force in the classical Darcy's law. The
second term accounts for flow due to a gradient in volume fraction. These terms
have previously been reported in [49]. The third is known as the Brinkman
correction factor and is often neglected for slow velocity flows. The coefficient
Ll is a fourthorder tensor resulting from the linearization process, and could be
a function of any independent variable that is not zero at equilibrium. The last
term on the first line accounts for flow due to thermal gradients and was also
reported in [49]; a comprehensive discussion of its impact can be found there.
The coefficient Hl is a thirdorder tensor arising from the linearization process
and could be a function of the same variables as Ll. It should be noted that
there are no isotropic thirdorder tensors, so that this term vanishes for isotropic
materials.
The first term on the second line accounts for gravitational effects. The
second term is previously unreported, although it resembles a term reported by
Singh et al. [49]. The firstorder component is responsible for flow induced by
the change in energy of the liquid phase due to (pure) shearing of the solid phase,
which is of interest at low moisture contents when only a few layers of water are
present. A similar term was originally reported by Murad and Cushman [40] and
Bennethum and Giorgi [11]. The secondorder and higherorder components are
(n)
novel. They are nonlinear in C and account for flow induced by the changes in
the liquid phase energy due to the rates of shearing. It is similar to dAl/d E'\
a term reported in [49], however Ea incorporates changes due to shear and
64
volumetric changes and it is unclear what volumetric changes are being captured
by Â£l and what by Es. Here we clearly distinguish between shearing effects and
changes due to moisture content, which are captured by the last term. The first
order component of the last term was previously reported by Achanta, [2], and
accounts for swelling in the normal directions and it has been speculated that
this term is responsible, in part, for nonFickian fluid transport. The second
and higherorder components of the last term are previously unreported, and
are nonlinear in time derivatives of volume fraction. As the polymer takes on
fluid the fluid acts as a plasticizer, lowering the glass transition temperature and
causing viscous relaxation and swelling. As these last terms are related to the
change in energy with respect to the time derivatives of the volume fraction of
the liquid phase, they are related to the normal components of strain; thus they
account for the effect of relaxation of the polymer matrix on fluid transport.
We now derive an equivalent form of equation (3.64) written in terms of the
Gibbs free energy which is related to the chemical potential of the species by
[13]
n i
Gl = =Al + ~v (3.65)
1=1 p
where the chemical potential is defined thermodynamically as [13]
d{sapAa)
d{eapai
= Aa + Pa
ea,T Op J
e",T
(3.66)
Writing (3.15) in terms of the independent variables plj = Cljpl (instead of
Clj and pl), we have
(")
s T=s .
A' = A'(Â£l, (eV>,r,C",C"),
(3.67)
65
where m 1,... ,p, n 1,...,q, and j = 1,... ,N. Thus, we have
s'p'VA1 = elp'
id A1
ds
__ / , ^
Ve' s'pl >  V e 1 eV VT
; i dA1
J=1
JL ("+
>6 JJ \ _LL_ (TT n
VplielplJ2
(>
=<9CS
+ 1)
: (V Cs )T.
(3.68)
Recall the definition of the swelling potential (3.43), and the three pressure
relationship (3.42). Combining these results to eliminate nl and then using this
to eliminate Â£lpl^r in (3.68) results in
JL f)Al ir i JL. 8A1 11111
fve Ey Â£ 4.vV' Â£y Â£ L:
m=ide>
'*= d Cs
N
pvA + p've' + Eyf vr + Â£>y~
j=l P
VpT (3.69)
The lefthand side of (3.69) can be substituted into equation (3.64) to obtain
8Al
Rl (eys) = elplVAl V(ey) + plVsl + elplVT
N
+Elpg + V (s'L' : d) + V (e'H1 VT) +
udA'
3 = 1
dp1^
Vp'T (3.70)
With the goal of writing equation (3.70) in terms of a gradient in the Gibbs free
energy, note that because pl} Cl] pl
Vp'J = plVCL + ClVp\
(3.71)
so that
N
Y,Â£lpl
3= 1
id A1
dp13
vpi>=YlÂ£l(p
i ]=\
i\2 dA^
dpli
VCl> + ^Vp',
P
(3.72)
where we have used the definition of the classical pressure, equation (3.35), with
a = l. Using (3.72) in (3.70), combining terms and rearranging, we obtain the
66
following form of Darcys Law
Rl (elv>s) = elplVGl + elpl^VT + Â£>V)2f4
j=i pJ
+Â£'plg + V(sl{Ll : d')) + VT)).
VC1'
(3.73)
The first term on the righthand side is the generalized Darcys law term which
says that flow is induced by a gradient in the Gibbs free energy. The second
term accounts for temperature effects due to rewriting Darcys law in terms of
the Gibbs potential (Gl is a function of temperature). The third term appears
because the gradient in the Gibbs free energy accounts for effects on the bulk
flow, but not entirely for contribution of species, which still may have an effect
on flow. The last three terms remain the same as in equation (3.64).
3.5 Transport Models
The goal is to derive a transport equation which can be used for swelling soils
and for swelling polymers. As such, we neglect higherorder terms and neglect
gravity since we expect the contribution of these terms to be small in comparison
with other terms, assume constant temperature, and assume negligible shearing
effects since in relevant experiments and applications fluid penetration occurs
in the normal directions. Thus, equations (3.64) and (3.73) simplify to the two
following equivalent forms of Darcys law
QAl ("+5
Rl (slvl's) = V(e'p') + plVsl elpl ^ (V Cs)
n=0 Q Cs
elpl
P dAl (mjj
VVT
m= ids1
R! (eV*) = eplVGl + Yzl{pl?^r
j=1 H
(3.74)
(3.75)
67
Either form of Darcys law, equation (3.74) or (3.75), can be used to derive a
transport equation. To avoid confusion we will call equation (3.74) the pressure
form of Darcys law, since the first term involves a pressure gradient, and we will
call equation (3.75) the Gibbs form. Both are valid for any twophase porous
medium that meets the requirements laid out in the introduction, where, as
mentioned above, gravitational and thermal effects have been neglected. In what
follows we restrict our derivation to materials in which each phase is composed
of only one species and where no mass transfer between phases takes place. For
such a system, the last term of equation (3.75) is zero. This means that the
lone term containing the gradient in Gibbs free energy accounts for all of the
information contained in the four terms of equation (3.74).
The same form of the continuity equations are used for both transport equa
tions. Assuming an incompressible liquid phase and rewriting the liquid phase
continuity equation in terms of Â£l Dssl/Dt we have
il + {l sl)V (sVs) = 0, (3.76)
where V vs was eliminated via V vs = il/es, which comes from equation (3.1)
assuming that the solid phase is incompressible.
3.5.1 Transport Equation in Terms of Pressures
The goal of this section is to obtain a transport equation that can be ex
pressed entirely in terms of the volume fraction of the liquid phase. To this end,
note that the first two terms on the righthand side of (3 74) can be written as
 V(eV) + f Vs' = el Vpl Â£l Vnl tt'Ve\ (3.77)
68
where we have used the three pressure relationship, (3.42). Now, we assuming
that nl is primarily a function of sl, that is
(3.78)
equation (3.77) becomes
V(eV) +piVe' = tt' + Â£lj elVpl.
(3.79)
Note that this assumption would not be valid for a nonswelling porous medium,
but it is reasonable for a porous medium composed of a highly interacting solid
and liquid phase.
Let K' be the unjacketed compressibility for a porous media with incom
pressible liquid and solid phases, [6]
where we have used equation (3.42) again. Substituting into (3.79) we obtain
(3.80)
p
then
(3.81)
(3.82)
where we have assumed pl = pl(sl, pl) to calculate Vp(. If we impose a constant
pressure along the boundaries then the density gradient Vp' will be negligible.
69
Neglecting Vpl and using (3.82) in the pressure form Darcys law, (3.74) we
have
slvls =
V
dAl
K'VeepÂ£ :(V C")
(n)
n=0 Q Cs
l _("
spy Vs:
m= 1 O Â£
(")(
,(3.83)
where we have assumed that the resistivity tensor Rl is isotropic, so that we
can write
Kl
(3.84)
where Kl is the permeability and p!v is the viscosity of the liquid phase. Substi
tuting elvl,s from (3.83) into (3.76) gives the first form of the general transport
equation
Â£l + (l
dAl,
K'Ve1 elpl V
= 0, (3.85)
m= 1 8 Â£ >
where we have chosen to neglect the terms involving C and its derivatives.
This means that the model will only be applicable to systems with relatively
slow flows with a moderate to high fluid content.
This closely resembles an equation obtained by Singh et al. in [51] for a three
scale system but with a few key differences. First, the equation obtained in [51]
contains a coefficient EwA, which is obtained in an ad hoc manner and is related
to the equilibrium elasticity of the material. Alternatively, in this formulation
we have used K', which is a macroscale physically measurable compressibility
possessing a rigorous thermodynamic definition. Second, the terms containing
appear naturally as a result of the constitutive independent variables;
there is no need to preform a change of variables from gradients in strain and
70
its time derivative to gradients in volume fraction and its time derivatives as in
(m)S
[51]. To make the change of variables from Es to Â£ Singh et al. assumed
that their system has moderate to high fluid contents. This assumption, as
mentioned before, means that shear forces are being neglected, and as this choice
of independent variables for the solid phase shows is equivalent to neglecting
terms involving Cf. Applying subsequent arguments of Singh et al. [51] to
(3.85), we obtain the following integrodifferential equation
El + (1 sl)V (^DVe1 Bv(t T)VÂ£ldT^j = 0, (3.86)
where
where
{Â£l)2KlK'
Â£sHlv
l2J
Bv[t) =
plKl
/A
m=1
(3.87)
(3.88)
d(m^S(t)
dt(rn~1)
(3.89)
is the time derivative of the dirac delta function, and where
dAl
^(m) o(m)(
O Â£ 1
(3.90)
are constant as long as A1 is a linear function of Â£l and its time derivatives. The
primary difference between (3.85) and the transport equation derived in [51] is
in the interpretation of coefficients D and Bv(t), and the ease with which this
equation is obtained. Appropriate choice of Bv(t r) allows this transport equa
tion to capture a wide range of viscoelastic material behavior. A comprehensive
discussion of its modeling capabilities can be found in [51].
71
3.5.2 Transport Model in Terms of the Chemical Potential
Next we turn to equation (3.75) and show how to combine it with the clas
sical FloryHuggins model used to model swelling polymers. To use equation
(3.75) we need to obtain an expression for the Gibbs free energy for the specific
system we wish to model. In this case we will seek an expression for the Gibbs
free energy that is valid for a polymer/solvent system. We use the standard
FloryHuggins theory for a polymer/solvent system, [26]. According to the the
ory the chemical potential for the solvent (water) and solute (polymer) phases
are given by
Hl = nl0 + RT[ln{ 1 us) + (1 l/x)va + x(^s)2] (3.91)
jis = fis0 + RT[ln(vs) + (x 1)(1 vs) + xz(l us)2], (3.92)
respectively, [26], Here, R is the gas constant and x is the ratio of molar volumes
of the solute and solvent. The interaction parameter, x, is as a measure of how
good a particular solvent is for a particular polymer, and is set to zero for good
solvents. The molar volume fraction parameters v1 and vs are given by
n
v =
nl + xns
xns
nl + xns
(3.93)
(3.94)
where nQ represents the number of aphase molecules. The units of vl are
von+vol Dividing the top and bottom by the total volume we find that vl
has the same units as sl. This suggests that we may take vl = sl. In fact, this
is precisely what is accomplished by using x, the ratio of molar volumes of the
solute to solvent in the denominator of (3.93) and (3.94), because x has the effect
72
of scaling the size of a solid phase molecule to take up the appropriate volume of
liquid phase molecules. Similarly, vs = es. Thus, we can write equations (3.91)
and (3.92) in terms of volume fraction
Hl = Mo + RT{ln{ 1 Â£s) + (1 l/x)Â£s + x(ss)2] (3.95)
// = /ig + RT[ln(ss) + (x 1)(1 es) + yx(l es)2]. (3.96)
Assuming a good solvent, we take X 0, and equations (3.95) and (3.96) become
fj! = fil0 + RT[ln(sl) + (1 \/x)es\ (3.97)
fis ns0 + RT[ln(ss) + (x l)sz]. (3.98)
For phases composed of only one species, the Gibbs potential of the phase is
equal to the chemical potential of the phase, so that we may take Gl = y} in
(3.75).
If x is sufficiently large, which is a valid assumption for polymer solvent
systems since the molar volume of most polymers is quite large in comparison
to the molar volume of solutes like water, equation (3.97) can be approximated
as follows
Hl Mo + ^[ln(s') + ss]
= fiq + RT[\n (el) + lneÂ£']
= Hl0 + RT[ In (s;e1_Â£,)l
= fil0 + RT[\n{a)}, (3.99)
where a Â£le1~Â£ is called the activity, and has been used in [4] to model a
penetrant into a glassy polymer. If the solution were instead composed of two
73
polymers we might take x ss 1. Substituting (3.99) into equation (3.75) we
obtain
Â£ivu = _^Tp'K\i _ (3.100)
Substituting equation (3.100) into (3.76) we obtain the following transport equa
tion
il + (1 el)V {DVs1) = 0,
(3.101)
where
D =
RTW
/4
i).
(3.102)
3.5.3 Model Comparison
The derivation of equations (3.86) and (3.101) used equivalent forms of
Darcys law and the same conservation of mass equations. However, the similar
ity between equations (3.86) and (3.101) is surprising because the development
of these equations came from two completely different perspectives. Equation
(3.101) simply lacks the integral term of (3.86), and the coefficient D of equation
(3.102) takes a slightly different form than that of equation (3.86). Moreover,
if we choose an appropriate form for Bv(t), as discussed in [51], equation (3.86)
reduces to equation (3.101) for polymers in the glassy state and in the rubbery
state. For polymers in the glassy or the rubbery state the decision to use the
coefficient associated with either equation should be determined by the applica
tion and the available experimental data for a given material. Equation (3.86)
is capable of capturing the complex behavior of polymers undergoing the transi
tion between the glassy state and rubbery state. The discussion of the necessary
experimental data to use equation (3.86) can be found in [53, 52].
74
3.6 Conclusions
In this chapter, we have shown that a judicious choice of independent vari
ables for the solid phase yields a physically meaningful expression for the solid
phase stress, validating the assumptions made by Bennethum and Weinstein in
[14]. We have shown how this choice of independent variables can be used to
elucidate the assumptions necessary to derive a transport equation applicable
to swelling viscoelastic materials analogous to that of Singh et al. [51], Using
FloryHuggins theory we derived another transport equation which is specifi
cally applicable to polymers. We then compared the resulting transport equa
tions (3.86) and (3.101) and argued that the latter should be thought of as a
simplification of the former to be used specifically when a polymer is in a glassy
or rubbery state, and the former should be used when the polymer is undergoing
a transition from one state to the other.
75
4. ThreePhase HMT for Swelling Porous Systems
Obtaining bulk and species transport equations that are capable of modeling
viscoelastic swelling systems in which mass transfer between phases takes place
requires special consideration. In this chapter we develop constitutive theory
for a system composed of a viscoelastic solid phase, an elastic solid phase, and
a viscous fluid that is capable of dealing with the unique difficulties that arise
in modeling such systems.
4.1 Motivation
In the context of HMT, fluid transport for swelling porous systems has been
investigated by a number of authors. Bennethum and Cushman presented theory
for such a system in [8]. In this work, they modeled the solid phase as elastic and
the liquid phase as viscous. The resulting model exhibited viscoelastic behavior
due to phase interactions. Singh et al. extended these results, modeling the
solid phase as viscoelastic, [49, 50]. These papers were aimed at biopolymeric
systems which display viscoelastic behavior at the microscale. In Chapter 3 we
revisited the transport theory for biopolymeric systems developed by Singh et
al. [51], in which they used the constitutive theory in [50] and the balance laws in
[7]. However, this transport theory did not include constituent transport, only
transport of bulk fluid, that is, transport of the liquid phase with respect to the
solid phase, but not transport of species within a phase. They also assumed
that no mass transfer takes place between phases. Because we wish to develop
76
a model that represents viscoelastic drug delivery systems, we must take into
account mass transfer between phases.
When attempting to model such systems using twophase HMT we en
counter a couple of difficulties. First, to obtain the bulk transport equation
(3.85) we assumed that both the solid and the liquid phases are incompressible,
that is
Unfortunately, because mass transfer between phases takes place, these condi
tions no longer hold. While we can make the argument that equation (4.2)
is approximately true by assuming that the dissolved drug in the liquid phase
is sufficiently dilute, we cannot make the same argument for the solid phase,
especially at high drug loadings.
Often, experiments in the literature measure the volume fraction of drug,
[20], In twophase HMT the variables that describe the constitution of the solid
phase are: the concentration of species within the phase, CS] = [mass s^/mass s],
the density of species within the phase, psi [mass Sj/vol s], the density of
the bulk phase, ps = [mass s/vol s], and the volume fraction of the phase,
Â£s [vol .s/vol]. Since none of these variables account for the volume of j in the
solid phase, twophase HMT lacks the ability to represent the volume fraction of
drug. This makes it difficult to compare the results of a twophase HMT model
to the available experimental data.
(4.1)
and
(4.2)
77
These considerations lead us to consider an alternative method of model
ing swelling drug delivery systems. In particular, a threephase model with a
viscoelastic solid phase (the polymer), an elastic solid phase (the drug), and a
viscous liquid phase, where we specify that mass transfer takes place between
the drug and liquid phases only, has the potential of overcoming the aforemen
tioned difficulties. This allows us to consider both solid phases as incompressible,
where equation (4.1) is the mathematical representation of this condition. We
can also argue that equation (4.2) is approximately true for the liquid phase, as
we did above. Furthermore, using three phases allows us to make more direct
comparisons with experiments since we now have a variable that represents the
volume fraction of drug in the solid phase.
In this chapter we will present the constitutive theory for such a system. The
details involved in deriving the entropy inequality are analogous to that found
in Chapter 3, so we will present only those details which differ significantly due
to the presence of three phases. The resulting entropy inequality will be given,
then we will present nonequilibrium, equilibrium, and nearequilibrium results,
as well as Darcys and Ficks laws.
4.2 Constitution
We assume the following independent variables:
l, Â£s, Vs1, Vss, (tl, V(tl, T, VT, pl, C
J /iL
()
Js% CCSv, C3v Cf\ JSv, Cs"i, va's, Va>, dl (4.3)
where m = l,... ,p denote material time derivatives of order m, n 1,...,q,
denote material time derivatives of order n, and j 1,...,JV 1 represent
78
N 1 constituents, l denotes the liquid phase, sv denotes the viscoelastic solid
phase, and se denotes the elastic solid phase. We do not include all three volume
fractions since they are related through
Â£l + Â£Sv + Â£Sr = 1, (4.4)
and therefore are not independent. We choose to use sl since it measures liquid
content of the system and this plays an important role in the swelling behavior
of the system. We choose to use e3v as opposed to Â£Sc because we want the
polymer phase to be the reference phase. We include T since the behavior of
most polymers is a function of temperature, and we include VT so that we
recover Fouriers law of heat conduction, that is, without it we obtain no heat
flux. The relative velocities va's and va^a give us information about exchange of
momentum within a phase and viscous diffusion. For the liquid phase we include
in our list of independent variables: pl, and Cl>, which account for the material
makeup, dl = 1 /2('Vvl + (Vvl)T) to account for the viscous nature of the liquid,
Vs1 to account flow due to moisture content, and Vr^1 to account for viscous
diffusion within the phase. For the viscoelastic solid we include: JSv and Csj
to account for the material makeup of the viscoelastic solid phase, derivatives of
order m = 1 to p for Â£l (derivatives in the volume fraction of the liquid phase are
related to the derivatives in the normal components of strain) and derivatives
of order 0 through q for CSv to account for (as well as JSv) the viscoelastic
nature of the material (combined they take the place of the usual higher order
derivatives in strain), Vs3 to account for flow due to polymer content, and
to account flow due to changes in moisture content. Finally, for the elastic
solid phase we include: JSe and CScj to account for material makeup, and Cs
79
(along with JSe) to account for elastic behavior.
4.3 Entropy Inequality
To simplify calculations, we postulate the dependence of the Helmholtz free
energies as follows:
(">
A1 = A1 (s', , s0*, pl, Cl* ,T,CS\ Cr),
(4.5)
()
As = As (Â£l, , (Â£);, J3, C3vi, T, cs, CSv ,cs"),
(4.6)
A3 = A3'(sl,e3v,lÂ£)l,J3',Csi
_ Jn)
T cs^cs
C
(4.7)
where Caj is the concentration of the jth species in phase a, and j = 1,..., N.
Otherwise we adhere to the axiom of equipresence: all other constitutive vari
ables are considered a function of list (4.3). Material time derivatives of A1, A*v,
and ASe appear in the entropy inequality (2.82). Using the chain rule, they may
be calculated as follows:
DlAl
Dt
dAl Dlsl dAl DlsSv
dsl Dt des" Dt
+ E
dA Dl{Â£h
g{,D Dt
dAl Dlpl
+ ~d^l)t
N1
8Al DlT dAl
+ 2_^phrv1 + +
i=i
Dt
dT Dt dC
DlCSv
Dt
+E
dAl
Ji*>
dcSv
(") ______________
Dl CSv dAl DlCSe
Dt h dCSe Dt
(4.8)
80
DS"AS' dAs" Ds*el dAs' Ds'es
Dt
dsl Dt
N1
+ E Fvi
+
dASv Ds"(Â£)l dASv Ds,,JSl1
E^^ +
m= 1
cte'*" Dt a'1; Dt dJs< Dt
m= 1 v/ t
Ds(Cs
j=i
Dt
+
(")
dT Dt
+
c>C
Dt
? &4 Â£) CT' &4S" DSC0
+E ; +
M
n=l Q q*>
Dt
dc
Dt
(4.9)
DSAS
Dt
dASc DSeel dA8Ds'es" dAs* D^t1 dA3DSe J3
+
del Dt deSv Dt
EU SA
N1
+ E^S
m= 1 <9T' Dt
DSc(CSej) dASc DSeT dASe D3cC
dJSc Dt
j=i
Dt
+
(n)
ar Dt
+
sc
Dt
+E
<9AS
(")
n=l Â£ c5*
Da C dA3' DSeC
+
Dt
dC
Dt
(4.10)
We can choose either solid phase to be the reference solid phase. We choose
the viscoelastic solid phase because in the case of drug delivery the elastic solid
phase will undergo phase transfer. As such, equation (3.23) is replaced by
DQ() DSv()
Dt
Dt
+ va^V()
(4.11)
Enforcing the conservation of mass equations using the Lagrange multipliers
given by equation (3.29) we arrive at the following form of the entropy inequality:
a=~4e ev
Q
1 (JJdAl
dAa \ DST
+ V
dT
Dt
f[Â£^ + sSvpS^ + sSep!
dA3 dA3e m
Se^SeplX j Â£l
81
dAs
* l Â£P'
, dAl
T \ de
N
+ Â£Slps
i d A*
des'
+ Â£Seps
des
Q j= 1
1 ( i /dA1 ;w \ /
f (p7wÂ£A I P
dp1
1_
~T
1
f
1_
T
dASr 1 gs,!
N
Â£ p
^psv
dJs 3 Js*
dASc 1 ga
Er(*
J=1
yv
d Js* 3 J
kEMt
j=i
Jfle
jsv
dAs
, ; <9A' s, s> 3ASc s s
P s 3 Â£ c p *' ;s f Â£ v p
dC dCe DC
N
fr^e n?e
Â£{fst1 {FST
1Se
: C
"f
, i dAl s 0A
gp _s 3 g cp'
<9C SC
dAs
Sc Ss+ Â£s' pSv=z
N
^(^T1 (E*^) (F^
2
i / N
dC
pSv s\Sv \&v __
~T Â£ P A . (FSv)T
vj = 1
: C
+f + :d'
4E(l>v"
a \ n 1
dAa \
()
9CS
_1
'r
g P
i ( 9 A1
dAl
Ws
+ V ) VT + gy^ : V(7Se +gy ( ^ ^ ) Vg'
dAl X1
8 Â£l Â£l
i ,dA s l t dAl <&, l(A 0AJ _(m)i
+e Ti^ Ve + ^ E ^ : V C +Â£ P E 7^ V e
n=0
dcv
=i d'ie 1
m= 1
TV
+^y (77 ~ 7)Vp/+^ (E plj xai)vch+K^
1 ^1
Se
V
1,3
82
1
T
Â£SepSe
^ vr + : VC8' + es'ps'^Vsl
DC
dsl
Â£S'pSc
+esys
dAs
CL 7?/lSe (") CL >9sc .
VC" + Â£ V' Â£ : V Cs" +< V Y1 V i
71 = 0
d C
i d Â£l
dAs' ASe A
+ ) vjse+Â£svse j^^ej XSejJ vScj
dpsc JSe J
+t;v + t;
v
1 p
?E
, ^A1 s s dAs s s dAs' \ (m+i,,
Â£ p + Â£Sc p e + Â£ p I Â£ 1
T ^ V ^ a<">' ()'
V Â£aVr f Q A r
+ T2 ] ^ + XT
l j = l
1
p<*3Vaa ( ^7 + ^(Va>'af ) t7 .
t7 ^tQ*   AaN) + /9q'(A XaN)
y2
a
a JV1
Q j= 1
+e ? e {Â£4r +~'> (*?'+n *"'* (Â£
a j=l v ^
+ (AQ AaK)V(Â£apa>) V [eV'(^Qjv Aa>)} l va*'a
VnQa
1 N
4e
r J
j=i
Aaj + ^Kja)2
Nl
4eee^
a j=1 /3^a
Aj + (uQ)2
2V
1
L
~y Xl/L'C
a /3^a
JV1
Aa + (vaSl )2 + AQ xaiC
3= 1
> 0.
(4.12)
where we have used identities (3.31) and (3.32) to remove the Nth component
dependence from terms involving vQj,a, and V^', [13], and used (3.26)(3.27)
to eliminate dSv and dSe, as we did in Chapter 3, in favor of the independent
83
variables C and Js, and C and Js,\ respectively.
4.3.1 Three Phase Pressures
Because our system now contains three phases instead of two, it is neces
sary to revisit the subject of defining pressures for the liquid and solid phases.
The changes in pressure definitions result from the fact that we now have two
independent volume fractions as opposed to one. For the liquid phase we make
the following definitions, s sv,se:
p(Â£,es,p!,C^ = (pl)2~
(4.13)
sl,es,C J
is the classical pressure. Because two out of three of the volume fractions are
held fixed and the sum of the volume fractions is one, that is, sl + eSe + sSv = 1,
we can write pl as a function of any two of the three volume fractions.
,dAl
= elp
i J
del
(4.14)
elpl3 ,es
is the thermodynamic pressure of the liquid phase holding the volume fraction
of the s phase fixed,
(4.15)
es ,pl ,cb
is the swelling potential of the liquid phase holding the volume fraction of the s
phase fixed. These pressures are related through:
pl(Â£l,es\pl,Cl)=pls(el1es,elpl>) + nls(el,es,pl,ClJ),
(4.16)
which is analogous to equation (3.42). In equations (4.5)(4.7) we chose the
Helmholtz potentials to be functions of el and eSv. Thus, we expect the thermo
dynamic pressure and swelling potential holding eSv fixed to appear. Another
84
important relationship between the swelling potentials exists
< eV
dAl
des
(4.17)
e'.p'.c'.J
This means that if the liquid phase and the se phase are noninteracting then
nlSv ~ tt1Sc, which implies = plSe because of equation (4.16).
For the solid phases, s = sv,se, we make the following definitions:
,dAs
ps(sl,es, Js, CSj) = psJs
dJs
(4.18)
Â£,es,Csi
as the classical pressure of the s phase. Analogous to p\ we note that p3 can be
defined using any two of the three volume fractions.
,<9AS
des
(4.19)
ea,J\Ci
is the thermodynamic pressure of the s phase holding the volume fraction of the
P phase fixed, where 8 / s.
.dAs
dJ3
(4.20)
c 3,p,C3i
is the swelling potential of the s phase holding the volume fraction of the (3
phase fixed, where fi ^ s. These pressures related through:
ps(sl,e3,J3, C3) = %{e0,e3,J3,C3>) + Tig(sa, Js,ps, C3j). (4.21)
which can be derived the same way as equation (3.46). We note that the differ
ence between equations (4.18) and (4.20) is in the variables that are being held
fixed. Again, we have a relationship between the swelling potentials
(4.22)
Â£s ,JS ,CSJ
ni na ~ s P
,dAs
ds13
85
where (3 ^ s, l, that is, /? is the other solid phase. This means that if the solid
phases are noninteracting then nf = irsg, which implies that pf = ps3 because of
equation (4.21).
Because the Helmholtz potential of the se phase is not a function of its own
volume fraction the pressure variable ps' does not appear naturally, and instead
we have terms involving (cb4s<,)/(<9Â£,)elv.. and (dAs')j{dÂ£Sv)\p Consider the
following two choices of independent variables for the se phase:
As'{el,Â£s\...) = ASe(Â£l,es,...). (4.23)
We note that because of the relationship between the volume fractions, equation
(4.4), both sets of variables contain the same information, ASe and As are
not two different potentials, but rather the same potential with two different
representations of the same set of independent variables. The total differential
of equation (4.23) is given by
DA8' =
dA8'
8e1
dASe
dÂ£ +
8e1
dAs
()8Sv
Se
j t dA
d>Â£ f
dÂ£Sc
dÂ£Sv + ...
dÂ£Se + . .
Then, for example
where
dA3* dASc dASe
0e1 0e1 esv dss' Â£se
(4.24)
(4.25)
dÂ£H"
dÂ£l
= l,
(4.26)
Â£,'SU
because of the relationship between volume fractions, equation (4.4). Similarly,
dAs
dÂ£s'
dASe
8eSc
(4.27)
86
Equations (4.25) allows us to make conversions such as
dsl
(4.28)
where we have used definition (4.19). In the following sections we will use
calculations such as these to simplify our results.
4.3.2 NonEquilibrium Results
The following variables are neither constitutive nor independent and appear
linearly in the entropy inequality:
where s = sVlse. This means that they can vary arbitrarily, and because these
variables appear linearly, to avoid violating entropy inequality (4.12), the coef
()+!)/
^ ,
(4.29)
ficients of these variables must be identically zero. This yields the following set
of nonequilibrium results:
(4.30)
(4.31)
(4.32)
(4.33)
(4.34)
87
(4.35)
N
j=l
a
(4.37)
(4.38)
Consider equation (4.30). A classical result for a single material was demon
strated in Chapter 2, and says that T and rj are dual variables. In this case we
will make the simplifying assumption that
Equation (4.31) gives an identity for the Lagrange multiplier and for the re
an identity for AQj. It has been used to obtain (4.35) and will also be used
throughout this chapter. Equations (4.33) and (4.34) are novel. They are a
direct consequence of the new choice of independent variables. Since it is rea
sonable to expect the diffusive velocities in the solid phases to be negligible, see
equation (A. 12), we can rewrite these two equations as:
Vq.
(4.39)
mainder of this section we will use it wherever X1 appears. Equation (4.32) gives
(4.40)
(4.41)
88
Equation (4.35) contains an Ns/lcomponent dependence which we will remove
later since we do not want a result which depends on the labeling of the con
stituents. Equation (4.36) is novel and allows us to solve for the Lagrange
multiplier ASc. Again assuming the diffusive velocities in the sephase are small,
we have
tSr = ptrX'I + 2(FSe) (FST (442)
Q
Taking 1/3 the trace of (4.42), using (4.40) to eliminate tse, and then solving
for ASe we obtain
A = H ^
ps 3 ^ Â£s* p3'dC
Â£apa dAa Se
7=Z7 C .
(4.43)
Substituting this back into equation (4.42) yields
= ~PSeI + K'e + ^Kk
2 v Â£apa 0Aa
3 Â£s<= dCSe '
(4.44)
where
f) A St f) 4 st,
til = 2ffF* ^ (Fsr + 2 (Fsf (4.45)
and
tI; = 2T" (4.46)
where tsa\ and t*seh are the Terzaghi and hydration stress tenors for the sephase,
respectively. Note that when we take one third the trace of (4.44) we recover
(4.40) so that the physical pressure in the se phase, l/3tr(ts<>), coincides with
the classical pressure of that phase, pSc.
89
In equation (4.37) and (4.38) Â£apa is begin held fixed, therefore we can bring
it inside the partial derivative. Doing so yields
dAf
and
()
ocSv
dAx
d el
= 0,
(4.47)
= 0,
(4.48)
where elplAl + Â£s^ps<ASc + Â£s'psASv = At is the total Helmholtz free energy.
This means that the total Helmholtz free energy is not a function of the qth
derivative of the right CauchyGreen tensor nor a function of the pth derivative
of the volume fraction of the liquid phase.
4.3.3 Equilibrium Results
We define that equilibrium to be when the following variables are zero:
il, is\ ~tS\ dl, bs vl's, vSc's\ {T\ VT, Vr'1, va>a, e^, eÂ£, (4.49)
where n = 1,.... q 1 and m = 1,... ,p 1. Using the method for obtaining
equilibrium results outlined in Chapter 2, we obtain the following results:
dAa
P =
dÂ£l
(4.50)
yfy^ = o,
' dsSv
(4.51)
tSv = psASv J + 2 F
v (sr^Â£aPa dAa \
(4.52)
90
