DIFFUSION OF A FLUID THROUGH A VISCOELASTIC SOLID

arXiv:1010.3488v2 [cs.CE] 8 Dec 2010

SATISH KARRA
Abstract. This paper is concerned with the diffusion of a fluid through a viscoelastic
solid undergoing large deformations. Using ideas from the classical theory of mixtures and a
thermodynamic framework based on the notion of maximization of the rate of entropy production, the constitutive relations for a mixture of a viscoelastic solid and a fluid (specifically
Newtonian fluid) are derived. By prescribing forms for the specific Helmholtz potential and
the rate of dissipation, we derive the relations for the partial stress in the solid, the partial
stress in the fluid, the interaction force between the solid and the fluid, and the evolution
equation of the natural configuration of the solid. We also use the assumption that the
volume of the mixture is equal to the sum of the volumes of the two constituents in their
natural state as a constraint. Results from the developed model are shown to be in good
agreement with the experimental data for the diffusion of various solvents through high
temperature polyimides that are used in the aircraft industry. The swelling of a viscoelastic
solid under the application of an external force is also studied.

1. Introduction
Several materials in the areas of of polymer mechanics, asphalt mechanics, and biomechanics that show non-linear viscoelastic behavior swell in the presence of a fluid. For
instance, polyimides which are known to show viscoelastic solid-like response (see Bhargava
(2007), Falcone and Ruggles-Wrenn (2009)) are used in the aerospace industry due to their
good performance at high temperatures (also see Ghosh and Mittal (1996)); these materials
in their service environment are known to swell in the presence of moisture. In addition,
asphalt based materials (that are well know to show non-linear viscoelastic fluid-like behavior) degrade in the presence of moisture (Kim et al. (2004)). Diffusion of biological fluids
through biological materials is another application wherein typically nutrition is provided
by the fluid that diffuses, and the amount of the stress or strain in the solid can control
the chemicals that are released (Rajagopal (2007)). Thus, there is a considerable interest
to understand how such viscoelastic materials deform and swell due to diffusion of a fluid.
Study of such a phenomenon is also of interest in geomechanics (Cohen (1992)) and food
industry (Singh et al. (2004)).
It is well known that the Darcy’s and Fick’s equations (Darcy (1856); Fick (1855)) that
are extensively used cannot predict swelling of the solid as well as the stresses in the solid.
In fact, Darcy’s equation is an approximation of the balance of linear momentum of the fluid
going through a rigid solid. To capture the swelling phenomena, several works have been
done using mixture theory (see review article by Rajagopal (2003)) and using variational
Date: March 17, 2018.
2000 Mathematics Subject Classification. 74D10, 74F20.
Key words and phrases. viscoelastic solid, mixture theory, maximization, entropy production, natural
configuration, polyimide.
Pre-print submitted to International Journal of Solids and Structures.
1

principles (Baek and Srinivasa (2004)). These models have been shown to match well with
experimental data for swelling of rubber (that show elastic response) due to the diffusion of
various organic solvents.
In the area of diffusion of a fluid through viscoelastic materials, some of the earliest works
were by Biot (1956) and Weitsman (1987), who have used linear viscoelasticity. In deriving
their models, they have used the fact that the fluxes and affinities are related through linear phenomenonological relations. Later on, Cohen and co-workers (Cohen and White Jr.
(1991); Edwards and Cohen (1995)), and Durning and co-workers (Huang and Durning (1997);
Cairncross and Durning (1996)) have also recognized the importance of studying the diffusion of solvents through polymeric materials that show viscoelastic response. They have
coupled diffusion and viscoelasticity by adding terms to the flux of the diffusing fluid that
depend on the stress in the solid. Such an approach does not have a thermodynamic basis. Furthermore, these models developed are only one-dimensional in nature. Recently,
Liu et al. (2005) have used the model developed by Cohen and co-workers to study the effect
of various viscoelastic parameters on diffusion. They have shown that comparable relaxation
times of polymer viscoelasticity and diffusion of a fluid results in non-Fickian behavior.
1.1. Main contributions of this work. Our main goal in this paper is to develop a
thermodynamic framework to model diffusion of a fluid through a viscoelastic solid and
we shall mainly focus on the swelling of polyimides. We use ideas from mixture theory
(see Truesdell et al. (2004), Truesdell (1957a), Truesdell (1957b), Atkin and Craine (1976),
Bedford and Drumheller (1983), Green and Naghdi (1969), Samohỳl (1987), Bowen (1982),
Rajagopal and Tao (1995) for details) and irreversible thermodynamics to build such a framework. In Karra and Rajagopal (2010), a framework that can be used to predict the non-linear
viscoelastic response of polyimides under various temperature and loading conditions, has
been developed. In the current paper, such a framework in Karra and Rajagopal (2010) is
extended to incorporate diffusion of a fluid and to model the swelling phenomenon.
The thermodynamic framework in the current work uses the notion of evolving natural configuration that has been used to model a variety of phenomena including classical plasticity, viscoelasticity, multi-network theory, superplasticity, twinning, etc. (see
Rajagopal and Srinivasa (2004a) for details). The evolution of such a natural configuration is
determined by maximizing the rate of entropy production (with any additional constraints).
We constitutively prescribe forms for the Helmholtz potential of the mixture and the rate of
dissipation (which is the product of density, temperature and the rate of entropy production)
due to mechanical working, diffusion and heat conduction. The final constitutive relations
are then derived by maximizing the rate of dissipation under appropriate constraints. In such
an approach, one need not assume linear phenomenonological relations between the flux and
the affinites, and thus our framework is more general. It has also been shown recently that
if one chooses a quadratic form for the rate of entropy production in terms of affinities, and
maximizes the rate of entropy production with respect to the affinities, one can arrive at the
Onsager’s relations (see Rajagopal and Srinivasa (2004b) for further details).
An initial boundary value problem is considered wherein a viscoelastic solid is held between
two rigid walls, and immersed in a fluid. Using the model developed in the current paper,
free swelling of such a solid and its swelling under the application of external force i.e., stressassisted swelling are studied. The numerical results for free swelling of the viscoelastic solid
2

are compared with experimental data for diffusion of different solvents through PMDA-ODA
and HFPE-II-52 polyimide resins.
1.2. Organization of the work. In section (2), the kinematics required in this paper are
documented. In section (3), the constitutive assumptions are specified and the constitutive
relations are derived. We shall also show that our constitutive relations reduces to the
equations for diffusion through an elastic solid derived using theory of mixtures when certain
parameters take special values. An initial boundary value problem is set up using our
model in section (4). The boundary conditions used, the non-dimensionalization scheme,
and comparison of the numerical results with experimental data are given in (4.1), (4.2),
and (4.3), respectively. Final concluding remarks are given in section (5).
2. Preliminaries
Let us consider a mixture of a Navier-Stokes fluid and a viscoelastic solid. We shall assume
co-occupancy of the constituents, which is the central idea in theory of mixtures. It is based
on the notion that at each point x in the mixture, at some time t, the two constituents exist
together in a homogenized fashion, and are capable of moving relative to each other. We
shall denote the quantities associated with the fluid through the superscript f and use the
superscript s for the solid. Now, we shall define the motion χi for the i-th constituent of the
mixture through

x = χi X i , t , i = f, s,
(2.1)

where X i is the material point of the i-th constituent in its reference configuration. We
shall assume that the mapping χi is sufficiently smooth and invertible at each time t. The
velocity associated with the i-th constituent is defined as
∂χi
vi =
,
(2.2)
∂t
and the deformation gradient through

∂χi
.
(2.3)
∂x
Let κt denote the current configuration of the mixture and let κR , κr denote the reference
configurations of the solid and the fluid respectively. Also, let κp(t) denote the natural
configuration of the viscoelastic solid (see Fig. (1)). Such a configuration is attained by the
solid body upon instantaneous removal of external loading. For a Navier-Stokes fluid, the
natural configuration is same as the current configuration of the fluid (Rajagopal (1995)).
Now, if F i , i = f, s is the gradient of the motion (usually known as deformation gradient)
χi X i , t , and if F sκp(t) is the gradient of the motion of the viscoelastic solid from κp(t) to
κt , then
−1

s
s
G = F κp(t)
F s.
(2.4)
Fi =

The density ρ and the average velocity (also known as barycentric velocity) v of the
mixture are defined by
X
1X i i
ρv.
(2.5)
ρi , v =
ρ=
ρ i
i
3

current configuration of the mixture

reference configuration of the solid

κR

κt

s

s

χ (X , t)
x
Xs

Gs

Fsκp(t)

κp(t)

χf (Xf , t)

natural configuration of the solid

κr
Xf

reference configuration of the fluid

Figure 1. Illustration of the various configurations of the viscoelastic solid
and fluid components in the mixture.
We define the following derivatives for any scalar quantity φi by
∂φi
∂ φ̃i
=
,
∂t
∂t

di φi
∂ φ̂i
=
,
dt
∂t

grad(φi ) =

∂ φ̃i
,
∂x

∂φi
∂ φ̂i
=
,
∂X i
∂X i

(2.6)

where
Hence,


φi = φ̃i (x, t) = φ̂i X i , t .

(2.7)

di φi
∂φi
=
+ grad(φi ) · v i ,
(2.8)
dt
∂t
and we shall also define the following
∂φ
dφ
=
+ grad(φ) · v.
(2.9)
dt
∂t
The velocity gradient for the i-th component Li and the velocity gradient for the total
mixture L are defined by
Li = grad(v i ),

L = grad(v).
4

(2.10)

The symmetric and anti-symmetric parts for the velocity gradients Li , L are
T i
T i
1h i
1h i
, Wi =
,
L + Li
L − Li
Di =
2
2
i
i
1h
1h
D=
L + (L)T , W =
L − (L)T .
(2.11)
2
2
The left Cauchy-Green stretch tensor B sp , B sp(t) and their principal invariants are defined as
T

s
s
s
s T
s
s
B G = F G (F G ) , B p(t) = F κp(t) F κp(t) ,
(2.12)

IBjs = tr B sj ,

IIBjs =

h
2
2 io
1 n
tr B sj
− tr B sj
,
2


IIIBjs = det B sj ,

j = G, p(t),
(2.13)

where det(.) is the determinant of a second order tensor. We shall next note the balance
laws.
The balance of mass for the i-th constituent without any mass production is given by

∂ρi
+ div ρi v i = 0,
(2.14)
∂t
where ρi is the mass density of the i-th constituent and div(.) := tr(grad(.)) is the divergence
operator with tr(.) meaning the trace of a second order tensor. The summation of Eq. (2.14)
over i along with Eq. (2.5) leads to
∂ρ
+ div (ρv) = 0.
(2.15)
∂t
The balance of linear momentum for i-th constituent is
h
T i
di v i
= div T i
+ ρi bi + mi ,
(2.16)
ρi
dt
where mi is the interaction force on the i-th constituent due to the other constituents, bi
is the external body force on the i-th constituent, T i is the partial Cauchy stress tensor
associated with the i-th constituent related to the surface traction on the i-th constituent ti
through
T
ti = T i n,
(2.17)
where n is the surface outward normal. From Newton’s third law, we have
X
mi = 0.

(2.18)

i

For mixtures, the balance of angular momentum, in the absence of body couples requires
that the total Cauchy stress of the mixture be symmetric i.e.,
X
T = T T , where T =
T i,
(2.19)
i

i

although the individual partial stresses T could be non-symmetric. Now, the balance of
energy for the i-th constituent is given by



v i .v i
i
i d
ǫ +
= div T i · v i − q i + ρi r i + ρi bi · v i + E i + mi · v i ,
(2.20)
ρ
dt
2
5

where ǫi , q i , r i are the specific internal energy, heat flux, radiant heating associated with
the i-th component and E i is the energy supplied to the i-th constituent from the other
constituents.
Now, taking the scalar multiplication of Eq. (2.16) with v i and subtracting the resulting
equation from Eq. (2.20), we arrive at
ρi

dǫi
= T i · Li − div(q i ) + ρi r i + E i ,
dt

(2.21)

Using ǫi = ψ i + θη i , where ψ i , η i are the Helmholtz potential and specific entropy of the i-th
constituent, and with θ being the common temperature of the constituents at a point in the
mixture, Eq. (2.21) along with Eq. (2.14) results in
 i
 1 i
∂ i i
1
q
1
1
i
i i i
ρ η + div ρ η v = T · L − div
− 2 q i · grad(θ) + ρi r i + E i
∂t
θ
θ
θ
θ
θ
 i i

i
i
dθ
ρ dψ
.
(2.22)
+ ηi
−
θ
dt
dt
i

Now, using the fact that η i = − ∂ψ
, we can establish the following result:
∂θ


 i

i
∂ψ i
di ψ i ∂ψ i di θ
∂ψ
∂ψ i ∂θ
di ψ i
i
i
id θ
+ v · grad(ψ ) −
+η
=
−
=
−
grad(θ)
dt
dt
dt
∂θ dt
∂t
∂θ ∂t
∂θ
 i i
dψ
,
(2.23)
=
dt θ fixed
where the subscript “θ fixed” means that the derivative is to be taken keeping θ fixed. We
shall define
X
1X i i
q=
qi, r =
ρr.
(2.24)
ρ i
i
Using the relation Eq. (2.23) in Eq. (2.22) and summing over i, along with Eq. (2.24), we
get
!
!
q
r
X
1X i i
∂ X i i
1
ρi η i v i =
T .L − div
− 2 q · grad(θ) + ρ
ρ η + div
∂t
θ i
θ
θ
θ
i
i


1 X i 1 X i di ψ i
+
E −
ρ
.
(2.25)
θ i
θ i
dt θ fixed
Eq. (2.25) is the balance of entropy with the rate of entropy production ζ being


1
1 X i 1 X i di ψ i
1X i
i
T · L − 2 q · grad(θ) +
E −
ρ
.
ζ=
θ i
θ
θ i
θ i
dt θ fixed

(2.26)

We shall assume that the total entropy production can be additively split into entropy production due to thermal effects i.e., conduction (ζc ), and entropy production due to internal
working and mixing (ζm ). We shall also require that each of these quantities be non-negative,
6

so that the rate of entropy production ζ is non-negative and the second law of thermodynamics is satisfied automatically. This implies that
q · grad(θ)
≥ 0,
θ2


1 X i 1 X i di ψ i
1X i
i
T ·L +
E −
ρ
≥ 0.
ζm :=
θ i
θ i
θ i
dt θ fixed
ζc := −

(2.27a)
(2.27b)

We shall choose q = −k(ρ, θ)grad(θ), k ≥ 0, so that Eq. (2.27b) automatically satisfies.
Also, if we define the rate of dissipation ξm := θζm , then
X
X
X  di ψ i 
i
i
i
T ·L +
E −
ρi
= ξm ≥ 0.
(2.28)
dt
θ
fixed
i
i
i
Assuming
X

Ei +

X

i

i

X

mi · v i = 0,

(2.29)

i

Eq. (2.28) can be re-written as
X

i

i

T ·L −

i

X

i

m ·v −

i

i

ρ

i



di ψ i
dt



= ξm .

(2.30)

θ fixed

Now,
X

ρi

i



di ψ i
dt



∂
=
∂t

X

ρi ψ i

i

!

X

+ div

ρi ψ i v i

i

!

!
X

dψ
=ρ
+ div
ρi ψ i v i − v ,
dt
i

P
where ψ := ρ1 i ρi ψ i is the average Helmholtz potential of the mixture.
Finally, from Eqs. (2.31) and (2.30), we arrive at
"
!#
X
X
X

dψ
T i · Li −
mi · v i − ρ
+ div
= ξm .
ρi ψ i v i − v
dt
i
i
i

(2.31)

(2.32)

θ fixed

Assuming that all the components have the same Helmholtz potential Eq. (2.32) reduces
to
X
i

i

i

T ·L −

X
i



dψ
m ·v − ρ
dt
i

i



= ξm ,

(2.33)

θ fixed

where we have used Eq. (2.5). The second law of thermodynamics is invoked by ensuring
ξm ≥ 0. The preliminaries discussed so far are sufficient for the derivation of the constitutive
equations in section (3).
7

3. Constitutive assumptions
We shall assume that the specific Helmholtz potential for the mixture is of the form


s
s
s
, IIBp(t)
, IIIBp(t)
ψ = ψ̂ θ, IBGs , IIBGs , IIIBGs , IBp(t)
,
(3.1)

and so

dψ
ds ψ
=
+ (v − v s ) · grad(ψ)
(3.2)
dt  dt
dψ
=
⇒
dt θ fixed
"
#
!
∂ ψ̂
∂
ψ̂
∂ ψ̂
∂ ψ̂
s
s
s
+ IBp(t)
I−
B sp(t) + IIIBp(t)
(B sp(t) )−1 · Ḃ p(t)
s
s
s
s
∂IBp(t)
∂IIBp(t)
∂IIBp(t)
∂IIIBp(t)
#
!
"
∂ ψ̂
∂ ψ̂
∂
ψ̂
∂ ψ̂
s
+ IBGs
I−
B sG + IIIBGs
(B sG )−1 · Ḃ G
+
∂IBGs
∂IIBGs
∂IIBGs
∂IIIBGs
+ (v − v s ) · (grad(ψ))θ fixed ,

(3.3)

ds ( )
where (˙ ) is
for the sake of convenience.
dt
Now,
s

s

s

Ḟ = Ḟ κp(t) Gs + F sκp(t) Ġ
s

s

s

⇒ Ḟ (F s )−1 = Ḟ κp(t) Gs (Gs )−1 (F sκp(t) )−1 + F sκp(t) Ġ (Gs )−1 (F sκp(t) )−1

In addition,

⇒ Ls = Lsp(t) + F sκp(t) LsG (F sκp(t) )−1
i
1h s
s
s
s
s T
s
s
s
−1
−T
T
F κp(t) LG (F κp(t) ) + (F κp(t) ) (LG ) (F κp(t) ) .
⇒ D = D p(t) +
2
s

s

(3.4)

s

Ḃ p(t) = Ḟ κp(t) (F s )T + F s (Ḟ κp(t) )T
= Lsp(t) B sp(t) + B sp(t) (Lsp(t) )T ,

(3.5)

and similarly
s

Ḃ G = LsG B sG + B sG (LsG )T .

(3.6)

Assuming that the response of the viscoelastic solid from the current configuration to its
natural configuration is isotropic elastic, we choose κp(t) such that
F sκp(t) = V sκp(t) ,
where V sκp(t) is the right stretch tensor in the polar decomposition of F sκp(t) .
8

(3.7)

Using Eqs. (3.4), (3.5), (3.6), (3.7) in Eq. (3.3), we get
 
dψ
=
dt θ fixed
"
!
#
∂ ψ̂
∂
ψ̂
∂
ψ̂
∂ ψ̂
s
s
2
B sp(t) −
(B sp(t) )2 + IIIBp(t)
+ IBp(t)
I · (Ls − LsG )
s
s
∂IsBp(t)
∂IIBp(t)
∂IIBp(t)
∂IIIsBp(t)
!
#
"
∂ ψ̂
∂
ψ̂
∂
ψ̂
∂ ψ̂
+ IBGs
(B s )2 + IIIBGs
B sG −
I · LsG
+2
∂IsBG
∂IIBGs
∂IIsBG G
∂IIIBGs
+ (v − v s ) · (grad(ψ))θ fixed ,

(3.8)

In what follows, we shall assume that the reference configurations (subscript o) of the
ρi
ρio
constituents are same as their natural states (subscript R) and so φi := i =
=
ρR
detF i ρiR
1
, i = s, f , where we have used the fact that ρio = ρiR . This need not be true in general.
detF i
We shall also assume the volume additivity constraint that is based on the fact that the
volume of the swollen viscoelastic solid is equal to the sum of the volumes of the unswollen
viscoelastic solid and the fluid (Mills (1966)). In our case this constraint is given by,
φs + φf = 1,

(3.9)

and so Eq. (2.14) can be re-written as

∂φi
+ div φi v i = 0,
∂t

which implies

!
i
X
φ
i
i
i
+ div
φv =0
∂t
i

⇒ div φf v f + φs v s = 0 (using Eq. (3.9)).

∂

P

(3.10)

(3.11)
(3.12)

Eq. (3.12) can be re-written as

φs tr(Ls ) + φf tr(Lf ) + v s · grad(φs ) + v f · grad(φf ) = 0.

(3.13)

Again from Eq. (3.9), we have
grad(φs ) + grad(φf ) = 0,

(3.14)

and hence, using Eq. (3.14) in Eq. (3.13), we arrive at
φs tr(Ls ) + φf tr(Lf ) + v s,f · grad(φs ) = 0,
where v s,f = v s − v f , is the velocity of the solid with respect to the fluid.
Next, we shall assume that the rate of entropy production is of the form

ξm = ξm LsG , F sp(t) , Lf , θ, v s,f ,

and so Eq. (2.33) along with Eq. (2.18) reduces to



dψ
s
s
f
f
s
= ξm LsG , F sp(t) , Lf , θ, v s,f .
T · L + T · L − m · v s,f − ρ
dt θ fixed
9

(3.15)

(3.16)

(3.17)

Using Eq. (3.8) in Eq. (3.17), we get
T s · Ls + T f · Lf − ms · v s,f − T sp(t) · (Ls − LsG ) − T sG · LsG − ρ(v − v s ) · (grad(ψ))θ fixed

= ξm LsG , F sp(t) , Lf , θ, v s,f ,
(3.18)

where

!
#
∂
ψ̂
∂
ψ̂
∂
ψ̂
∂
ψ̂
s
s
T sp(t) := 2ρ
+ IBp(t)
B sp(t) −
(B sp(t) )2 + IIIBp(t)
I , (3.19)
s
s
s
∂IBp(t)
∂IIBp(t)
∂IIBp(t)
∂IIIsBp(t)
"
!
#
∂
ψ̂
∂
ψ̂
∂
ψ̂
∂
ψ̂
T sG := 2ρ
+ IBGs
B sG −
(B sG )2 + IIIBGs
I .
(3.20)
∂IBGs
∂IIBGs
∂IIBGs
∂IIIBGs
"

Eq. (3.18) with the constraint Eq. (3.15) can be written as
T s · Ls + T f · Lf − ms · v s,f − T sp(t) · (Ls − LsG ) − T sG · LsG − ρ(v − v s ) · (gradψ)θ fixed


+ λ φs tr(Ls ) + φf tr(Lf ) + v s,f · grad(φs ) = ξm LsG , F sp(t) , Lf , θ, v s,f ,
(3.21)

where λ is a Lagrange multiplier.
We shall further assume that the rate of dissipation can be additively split into the rate of
dissipation due to mechanical working of the viscoelastic solid, the rate of dissipation due to
the Navier-Stokes fluid and the rate of dissipation due to diffusion of the fluid, with specific
forms as follows:


ξm LsG , F sp(t) , Lf , θ, v s,f = ξ LsG , B sp(t) , θ + νD f · Df + α(θ)v s,f · v s,f .
(3.22)
Then, from Eq. (3.22) and Eq. (3.21), we arrive at

T s · Ls + T f · Lf − ms · v s,f − T sp(t) · (Ls − LsG ) − T sG · LsG − ρ(v − v s ) · (grad(ψ))θ fixed


+ λ φs tr(Ls ) + φf tr(Lf ) + v s,f · grad(φs ) = ξ LsG , B sp(t) , θ + νD f · Df + α(θ)v s,f · v s,f ,
(3.23)
which can be re-written as





Ls · T s + λφs I − T sp(t) + Lf · T f + λφf I − νD f + T sp(t) − T sG · LsG



+ v s,f · −ms + λ grad(φs ) − α(θ)v s,f + ρf (grad(ψ))θ fixed = ξ LsG , B sp(t) , θ ,

(3.24)

using the fact that ρ (v − v s ) = −ρf v s,f . Since, the right hand side of Eq. (3.24) does not
depend on Ls , Lf and v s,f , we have
T s = −λφs I + T sp(t) ,
f

(3.25)

f

f

T = −λφ I + νD ,
s

(3.26)

s

f

m = λ grad(φ ) − α(θ)vs,f + ρ (grad(ψ))θ fixed ,

(3.27)

and so Eq. (3.24) reduces to


T sp(t) − T sG · LsG = ξ LsG , B sp(t) , θ .

(3.28)

Next, we shall maximize the rate of dissipation ξ with Eq. (3.28). We shall maximize the
auxiliary function



Φ := ξ + β ξ − T sp(t) − T sG · LsG .
(3.29)
10

Now,

∂Φ
(β + 1) ∂ξ
s
s
s = 0 ⇒
s − T p(t) − T G = 0.
∂LG
β ∂LG

(3.30)

Taking the scalar product of Eq. (3.30) with LsG and using Eq. (3.28), we arrive at
(β + 1)
ξ
=
,
∂ξ
β
s
· LG
∂LsG

(3.31)

and hence the evolution equation for the natural configuration of the solid is given by

T sp(t) − T sG =

ξ
∂ξ
· LsG
∂LsG

∂ξ
.
∂LsG

(3.32)

3.1. Specific constitutive assumptions. We shall assume the following specific form for
the specific Helmholtz potential of the mixture
 
µG0 − µG1 θ
cs1
θ
2
s
s
s
s
+
(IBGs − 3)
ψ̂ = A + (B + c2 ) (θ − θs ) − (θ − θs ) − c2 θln
2
θs
ρs θs

Rθ 
µp0 − µp1 θ
s
s
s 2
s
(1
−
φ
)ln(1
−
φ
)
−
χ(φ
)
,
(3.33)
(I
−
3)
+
+
B
p(t)
ρs θs
ρfR V0 φs

where µG0 , µG1 , µp0, µp1 are material parameters, θs is a reference temperature for the viscoelastic solid, R is the gas constant, θ the absolute temperature of mixture, V0 the molar
volume of the fluid, and χ a mixing parameter for the particular solid-fluid combination. The
last term in Eq. (3.33) is the term added to the specific Helmholtz in (Karra and Rajagopal
(2010)) to capture the swelling phenomenon in the solid.
Now,
η=−

∂ ψ̂
∂θ
s

= −(B +
−

cs2 )

R
ρfR V0 φs



+

cs1

(θ − θs ) +

cs2 ln



θ
θs



+ cs2 +

µG1
µp1
(IBGs − 3) +
(IBs − 3)
ρθs
ρθs p(t)


(1 − φs )ln(1 − φs ) − χ(φs )2 .

(3.34)

The internal energy ǫ is given by
ǫ = ψ + θη

= As − B s θs + cs2 (θ − θs ) +

 µG0
µp0
cs1 2
θ − θs2 +
(IBGs − 3) +
(IBs − 3),
2
ρθs
ρθs p(t)

(3.35)

and the specific heat capacity Cv is

Cv =

∂ǫ
= cs1 θ + cs2 .
∂θ
11

(3.36)

From Eq. (3.33) and Eq. (3.19)


 i
ρJps JGs h
s
s
s
s
T p(t) =
2µ̄p B p(t) + µ̄p IBp(t) − 3 I + µ̄G IBG − 3 I
ρsR

ρRθJps JGs 
s
s
s 2
ln(1
−
φ
)
+
φ
+
χ(φ
)
I,
+
ρfR V0

and from Eq. (3.33) and Eq. (3.20),


 i
ρJps JGs h
s
s
s
T sG =
2µ̄
B
+
µ̄
I
−
3
I
+
µ̄
−
3
I
I
G G
p
Bp(t)
G
BG
ρsR

ρRθJps JGs 
s
s
s 2
ln(1
−
φ
)
+
φ
+
χ(φ
)
I,
+
ρfR V0
µ̄

(3.37)

(3.38)

θ

G1 θ
, µ̄p = p0 θs p1 , JGs = det(Gs ), Jps = det(F sκp(t) ).
where µ̄G = µ̄G0 −µ̄
θs
We shall further assume that the rate of dissipation ξ is of the form

−µ̄

ξ = γ(θ)D sG · DsG ,

(3.39)

then Eq. (3.25) becomes


 i
ρ h
s
s
2µ̄
B
+
µ̄
I
−
3
I
+
µ̄
−
3
I
I
p p(t)
p
Bp(t)
G
BG
ρs

ρRθJps JGs 
ln(1 − φs ) + φs + χ(φs )2 I,
+
f
ρR V0

T s = −λφs I +

(3.40)

and Eq. (3.32) reduces to


2ρ 
s
s
µ̄
B
−
µ̄
B
= γ(θ)D sG .
p
G
p(t)
G
ρs

The final constitutive equations are


 i
ρ h
s
T s = −λφs I + s 2µ̄p B sp(t) + µ̄p IBp(t)
− 3 I + µ̄G IBGs − 3 I
ρ
s s 

ρRθJp JG
+
ln(1 − φs ) + φs + χ(φs )2 I,
f
ρR V0

(3.41)

(3.42a)

T f = −λφf I + νD f ,

(3.42b)

ms = λ gradφs − α(θ)v s,f + ρf (gradψ)θ fixed ,

(3.42c)

and

2ρ 
µ̄p B sp(t) − µ̄G B sG = γ(θ)D sG ,
s
ρ

(3.43)

being the evolution equation of the natural configuration of the solid.
Notice that when µ̄G = 0 and γ → ∞, for the LHS of Eq. (3.43) to be finite, we must have
D sG = 0. This implies that Gs = I, and hence B sp(t) = B s and the solid is now an elastic
12

solid. In such a case, the constitutive equations Eq. (3.42), with additional assumption of
ν = 0, reduce to
"
!
#
∂
ψ̂
∂
ψ̂
∂
ψ̂
∂
ψ̂
T s = −λφs I + 2ρ
Bs −
+ IB s
(B s )2 + IIIBs
I ,
(3.44a)
∂IBs
∂IIBs
∂IIBs
∂IIIBs
T f = −λφf I,

(3.44b)

ms = λ gradφs − α(θ)v s,f + ρf (gradψ)θ fixed .

(3.44c)

These equations are same as the equations derived using theory of mixtures for the diffusion
of a fluid through an elastic solid (see equations 3.15–3.17 in (Hron et al. (2002))).
4. Initial boundary value problem
Let us consider the problem of compression of the viscoelastic solid body inside rigid walls
as shown in Fig. (2). Let us assume that the motion of the swollen solid body is given by
x = X,

y = Y,

z = f (Z, t).

(4.1)

In this case, the deformation gradient of the solid (F s ) is given by
F s = diag {1, 1, p} ,
where p :=

(4.2)

∂f
. Let us assume a form for Gs as follows
∂Z
Gs = diag {1, 1, g(Z, t)} ,

(4.3)

and the velocity of the fluid be of the form
v f = (0, 0, v(Z, t)) .

(4.4)

(

(4.5)

Then,
B sκp(t)

s
IBp(t)

 2 )
p
= diag 1, 1,
,
g

 2
p
=2+
,
g

IBGs = 2 + g 2 ,

p
Jps = ,
g

JGs = g.

(4.6)

The balance of mass for the solid gives
φs =

1
1
s = ,
det(F )
p

(4.7)

and hence from volume additivity constraint
f

ρ =

ρfR



1
.
1−
p
13

(4.8)

Also, we note the following relations
ρ
ρfR
=
1
+
(p − 1),
ρs
ρsR
ρ
1
=1+
f
p
ρR

(4.9a)

!
ρsR
−1 .
ρfR

(4.9b)

The balance of mass for the fluid reduces to

fluid

F (t)

Porous filter

Z=H
Swollen viscoelastic solid
Z=0

Z = −H
F (t)
Figure 2. Schematic of the initial boundary value problem

∂v
∂p v ∂p
+
+ (p − 1)
= 0.
∂t p ∂Z
∂Z

(4.10)

Setting ν = 0, zz-component of the total stress tensor for the mixture reduces to
s
f
Tzz = Tzz
+ Tzz

(4.11)
!"

 2
 2


p
p
ρfR
2
2µ̄p
+ µ̄p
−
1
+
µ̄
g
−
1
= −λ + 1 + s (p − 1)
G
ρR
g
g2
!!
 


1 ρsR
1
1
Rθp
1
+ 1+
−1
ln 1 −
+ +χ 2 .
p ρfR
V0
p
p
p

#
(4.12)

The balance of linear momentum for the solid and the fluid after assuming zero body forces
reduce to
s
∂2f
∂Tzz
=
+ msz ,
∂t2
∂z
f
∂v
∂Tzz
ρf
=
− msz .
∂t
∂z

ρs

14

(4.13a)
(4.13b)

Now, adding Eqs. (4.13a), (4.13b), we arrive at the balance of linear momentum for the
mixture
∂2f
∂v
∂
f
ρs s + ρf
=
(T s + Tzz
)
∂t
∂t
∂z zz
2
sf
∂λ ∂Tzz
s∂ f
f ∂v
⇒ρ
+ρ
=−
+
,
(4.14)
∂ts
∂t
∂z
∂z
where
#
!"

 2
 2
f

ρ
p
p
sf
Tzz
= 1+ R
− 1 + µ̄G g 2 − 1
(p − 1)
2µ̄p
+ µ̄p
ρsR
g
g2
!!
 


1 ρsR
1
1
1
Rθp
+ 1+
ln 1 −
+ +χ 2 .
−1
(4.15)
p ρfR
V0
p
p
p
Now, Eq. (4.13b) along with volume additivity constraint reduces to
 


∂ψ
∂f
f ∂λ
f
f ∂v
,
= −φ
+α
−v −ρ
ρ
∂t
∂z
∂t
∂z θ fixed

(4.16)

∂f
).
∂t
Multiplying Eq. (4.14) with φf and subtracting Eq. (4.16) from the resulting equation we
get
 


2
sf
∂ψ
∂f
f s∂ f
f
f ∂v
f ∂Tzz
f
,
(4.17)
φ ρ
+
(φ
−
1)ρ
=
−α
−
v
+
φ
+
ρ
∂t2
∂t
∂t
∂z
∂z θ fixed
where we have used the fact that the velocity of the solid is (0, 0,

which reduces to
2

∂v
f
+ (φf − 1)ρf
= −α
φ ρ
2
∂t
∂t
f s∂




∂T sf
∂f
− v + φf zz + ρf
∂t
∂z

∂ ψ̃
∂z

!

,

(4.18)

where
 pµ̄p
µ̄G p
ψ̃ = s g 2 − 1 + s
ρR
ρR




 



1
1
1
Rθp
p2
1−
ln 1 −
−χ 2 .
−1 + f
g2
p
p
p
ρR V0

(4.19)

Now, assuming that the velocity and acceleration of the solid are small compared to that of
∂2f
∂f
and 2 in Eq. (4.18), we get
the fluid, we shall drop
∂t
∂t
!
sf
∂T
∂
ψ̃
∂v
zz
.
(4.20)
= αv + φf
+ ρf
(φf − 1)ρf
∂t
∂z
∂z
Next, we shall also assume that the acceleration of the fluid is also small and we shall drop
∂v
term in Eq. (4.20), to get
∂t
"
!#
sf
1 f ∂Tzz
∂
ψ̃
φ
.
(4.21)
v=−
+ ρf
α
∂z
∂z
15

Using Eq. (4.21) in Eq. (4.10), we arrive at
"
!#
"
sf
sf
∂p
∂T
p
−
1
1 ∂p
ρf
∂
ψ̃
∂
φf ∂Tzz
zz
f
f
φ
+
=
+
ρ
+
∂t
αp2 ∂Z
∂Z
∂Z
α ∂Z p ∂Z
p

∂ ψ̃
∂Z

!#

.

(4.22)

Also, note that


1 ∂g
D G = LG = diag 0, 0,
g ∂t



,

and so the evolution equation of the natural configuration reduces to
#
!"  
2
f
1 ∂g
ρ
p
γ
− µ̄G g 2 .
= 2 1 + sR (p − 1)
µ̄p
g ∂t
ρR
g

(4.23)

(4.24)

4.1. Boundary conditions. Applying boundary conditions in an initial boundary value
problem has been an issue in mixture theory due to its basic assumption of co-occupancy. For
instance, if traction is applied on the boundary, a natural question is how is the traction to be
split between the solid and the fluid. To this end, the method of spitting the traction based
on the volume fraction of the solid and the fluid was proposed (Rajagopal and Tao (1995)).
Later on, Baek and Srinivasa (2004) derived the relations for swelling of an elastic body based
on variational principles and the boundary conditions were derived naturally. However, this
approach assumes that the swelling is slow, and that the relative velocity between the solid
and the diffusing fluid is small. Recently, Prasad and Rajagopal (2006) have compared the
solutions of diffusion of a fluid through a elastic slab using various boundary conditions like
saturation boundary condition, traction splitting boundary condition, the natural boundary
condition derived by Baek and Srinivasa, and the condition that the chemical potential is
continuous across the boundary. Interestingly, they show that the results are insensitive to
these different forms of boundary conditions.
For our problem, let F (t) be the compressive force applied on the solid at Z = ±H as
shown in Fig. (2) and let P∞ be the pressure in the fluid at the boundaries Z = ±H, then
we shall apply the following boundary conditions:
s
Tzz
= −F (t) − φs P∞ ,
f
Tzz
= −φf P∞ ,

Z = ±H,

Z = ±H,

(4.25a)
(4.25b)

that is, we are assuming that the external force is borne by the solid only, while the fluid
pressure is borne by both the solid and the fluid, and this pressure is split proportional to
the volume fraction of the constituents. Based on these assumptions, Eq. (4.25b) reduces to
λ = P∞ ,

Z = ±H.

(4.26)

Eq. (4.25a) and Eq. (4.26) reduce to
sf
−F (t) = Tzz
,

Note that F (t) is zero under free-swelling.
16

Z = ±H.

(4.27)

4.2. Non-dimensionalization. We shall use the following non-dimensionalization scheme:
Z⋆ =

Z ⋆
t
vT ⋆
µ̄p
γV
, t = , v⋆ =
, p = p, g ⋆ = g, µ̄⋆p = , γ ⋆ =
,
L
T
L
µ
µL

where T , L are characteristic time and length respectively. If we pick µ =

γ⋆

∂p⋆
∂t

1 ∂g ⋆
g ⋆ ∂t⋆

L2 V0 α
RθT

Rθ
V0

and define the

, β2 :=
then Eqs. (4.22), (4.24) become
"
#


sf ⋆
1 ∂p⋆ ∂Tzz
1
∂ ψ˜⋆
= ⋆ 2 1− ⋆
+
(p )
p
∂Z ⋆ ∂Z ⋆
∂Z ⋆
"
!#

sf ⋆
˜⋆
∂
1
∂
ψ
1
∂T
zz
+ (p⋆ − 1) ⋆
+
,
(4.29a)
1− ⋆
∂Z
p
p⋆
∂Z ⋆
∂Z ⋆
#
 "  ⋆ 2

p
1 ⋆
= 2 1 + (p − 1) µ̄⋆p
(4.29b)
− µ̄⋆G (g ⋆)2 ,
β1
g⋆

non-dimensionalization quantities β1 :=
β2

ρsR
ρfR

(4.28)

where
#
 ⋆ 2
"

 ⋆ 2


p
(p )
1 ⋆
⋆
⋆
⋆ 2
⋆
− 1 + µ̄G (g ) − 1
+ µ̄p
= 1 + (p − 1) 2µ̄p
β1
g⋆
(g ⋆ )2

  


1
1
1
1
⋆
+ 1 + ⋆ (β1 − 1) p ln 1 − ⋆ + ⋆ + χ ⋆ 2 ,
(4.30a)
p
p
p
(p )


 



 p⋆ µ̄⋆p (p⋆ )2
1
1
1
µ̄⋆G p⋆
⋆ 2
⋆
⋆
˜
1 − ⋆ ln 1 − ⋆ − χ ⋆ 2 .
(g ) − 1 +
−1 +p
ψ =
β1
β1
(g ⋆ )2
p
p
(p )
(4.30b)

sf ⋆
Tzz

The boundary conditions Eq. (4.27) reduce to
sf ⋆
−F ⋆ (t⋆ ) = Tzz
,

Z ⋆ = ±H ⋆ ,

(4.31)

⋆

⋆
with F ⋆ = Fµ , P∞
= Pµ∞ . If we pick H as the characteristic length L, then H ⋆ = 1. The
coupled equations Eqs. (4.29a), (4.29b) are solved using the staggered scheme shown in
algorithm 1.
The ratio of the mass of the swollen solid to its original unswollen mass can be calculated
as follows:
R
ρ dV
m
=R s
m0
ρR dV
R z=1 ρ
dz
z=−1 s
ρR
= R z=1
dz
z=−1

Z z=1 
1
1
1 + ⋆ (β1 − 1) dz.
(4.32)
=
2β1 z=−1
p

Once, the value of p⋆ (Z ⋆ , t⋆ ) is evaluated on the domain at various times, Eq. (4.32) is
integrated numerically to get the mass ratio.
17

Algorithm 1 A staggered procedure for solving the coupled equations
⋆
1: Input: β1 , β2 , χ, F ⋆ , P∞
, µ̄⋆p , µ̄⋆G , γ ⋆ ; Time of integration t⋆f ; Time step ∆t⋆ ; No. of
⋆
divisions along Z direction N; TOLERANCE.
2: Output: p⋆ .
3: Set p⋆0 = 1, g ⋆0 = 1.
4: while t < tf do
5:
t⋆ = t⋆ + ∆t⋆ .
6:
while true do
⋆
7:
Using p⋆i(l) , g ⋆i(l) , N, ∆t, β1 , β2 , χ, F ⋆ , P∞
, µ̄⋆p , µ̄⋆G , γ ⋆ , Solve for p⋆i(l+1) using
Eq. (4.29a) in MATLAB’s pdepe solver with the p⋆i(l+1) at the boundaries (Z ⋆ = ±1)
obtained by solving the non-linear algebraic equation in Eq. (4.31). This non-linear
algebraic equation is solved using the fsolve solver in MATLAB.
⋆
8:
Using p⋆i(l+1) , g⋆i(l) and ∆t, β1 , β2 , χ, F ⋆ , P∞
, µ̄⋆p , µ̄⋆G , γ ⋆ , Solve for g ⋆i(l) using
Eq. (4.29b) in MATLAB’s ode45 solver.
9:
if kp⋆i(l+1) − p⋆i(l) k2 < TOLERANCE then
10:
Return.
11:
end if
12:
end while
13:
p⋆i+1 ← p⋆i , g ⋆i+1 ← g ⋆i .
14: end while
4.3. Comparison with experimental data. Fig. (3) shows comparison of the numerical results to the experimental data for the ratio of swollen to unswollen (see Eq. (4.32))
PMDA-ODA (poly(N, N’- bisphenoxyphenylpyromellitimide)) due to diffusion of the solvents DMSO (dimethylsulfoxide) and NMP (N-methyl-2-pyrollidinone). In case of DMSO
diffusing through PMDA-ODA the following material parameters were chosen: Density of
DMSO was chosen to be 1.096 g/cc (Yang et al. (2008)) and density of PMDA-ODA to be
1.42 g/cc (Srinivasan et al. (1994)), and so β1 = 1.3. Also, χ = 0.425, β2 = 0.018, µ̄⋆p = 0.1,
µ̄⋆G = 0.1, γ ⋆ = 20 were chosen. The characteristic time (T ) chosen was 10500 min. For
the diffusion of NMP, the material parameters chosen were: density of NMP is taken to
be 1.02 g/cc (Yang et al. (2008)) and so β1 = 1.4. Next, χ = 0.6, β2 = 0.016, µ̄⋆p = 0.1,
µ̄⋆G = 0.1, γ ⋆ = 20 and the characteristic time chosen was 245 min. The numerical results
show good agreement with the experimental data taken from (Gattiglia and Russell (1989)).
Next, we shall consider the diffusion of water through HPFE-II-52. The material parameters were assumed to be: χ = 0.425, β1 = 1.3, β2 = 0.018, µ̄⋆p = 0.1, µ̄⋆G = 0.1, γ ⋆ = 20. The
characteristic time chosen was 2800 s. Even in this case the numerical results and experimental data taken from (Bhargava (2007)) match well (see Fig. (4)). In all the numerical
calculations TOLERANCE was chosen to be 10−4 .
We shall now consider the problem of compression of the viscoelastic solid and study its
effects on swelling due to diffusion of a fluid. In this numerical experiment, the solid is
allowed to swell freely first till it saturates with fluid (upto t⋆ = 0.5). Then, the swollen solid
is subjected to constant compressive force of F ⋆ = 1 is applied for a time period of t⋆ = 0.5
and then the load is removed, and the solid is allowed to swell freely again for another time
period of t⋆ = 0.5. Fig. (5a) shows that the volume of the solid gradually increases with time
and then reaches a steady state where the volume of the solid is same everywhere and there
18

1.35
1.3
1.25

m
m0

1.2
1.15
1.1
1.05
1

Model
Experiment
0

1000

2000

3000

4000

5000

t(min)
(a) DMSO
1.25

1.2

m
m0

1.15

1.1

1.05
Model
Experiment
1

0

20

40

60

80

100

t(min)
(b) NMP

Figure 3. Comparison of the model with the experimental data from
(Gattiglia and Russell (1989)) for the diffusion of DMSO and NMP through
PMDA-ODA (imidized at 300oC) under free-swelling condition. The characteristic times chosen were 10500 min and 245 min for DMSO and NMP, respectively. Here, 301 spatial points were used for the calculations, non-dimensional
time step chosen is ∆t⋆ = 0.025.

is no further swelling. Also, the volume of the solid near the boundary increases faster than
that of the inner solid. Upon application of a constant compressive load in Fig. (5b), the
fluid diffuses out of the swollen solid and the volume of the solid gradually decreases until
the volume of the solid is same everywhere. Next, upon removal of the compressive load in
19

Normalized mass

1
0.8
0.6

0.4
0.2
Model
Experiment
0

0

200

400

600

800

1000

1200

1400

t(s)

Figure 4. Comparison of the model with the experimental data from
(Bhargava (2007)) (pg. 27) for the diffusion of water through HFPE-II-52
under free-swelling condition. The characteristic time chosen was 2800 s. The
parameters chosen are χ = 0.425, β1 = 1.3, β2 = 0.018, µ̄⋆p = 0.1, µ̄⋆G = 0.1,
γ ⋆ = 20. Here, 301 spatial points were used for the calculations, ∆t⋆ = 0.025.
m(t) − m0
, where m0 is the mass of the
The normalized mass is defined by
m∞ − m0
dry solid, m∞ is the steady state mass of the swollen solid, m(t) is the mass
of the swollen solid at a given time t.
Fig. (6), the solid absorbs the fluid and swells freely back to its original swollen saturation
state.
5. Conclusions
We developed a systematic framework with a thermodynamic basis to develop constitutive
relations for the diffusion of a fluid through a viscoelastic solid. A model was also derived
using this framework by choosing specific forms for the Helmholtz potential and the rate of
dissipation, and by maximizing the rate of dissipation. An initial boundary value problem
was solved where we considered free swelling and swelling under the application of external
force. We also showed that the model fits well with the experimental data for free swelling
of polyimides. Furthermore, our work in this paper can be easily extended to study the
diffusion of fluids in biological materials as well as in studying moisture-induced damage in
asphalt mixes and other geomaterials that show viscoelastic behavior.
Finally, here we shall summarize the assumptions made in this paper:
(1) the specific Helmholtz potential of the constituents is the same,
(2) the temperature of the constituents is the same,
(3) the specific Helmholtz potential of the mixture depends on the temperature of the
mixture, and the deformation of the solid,
20

1.4
1.35
1.3
1.25

J

1.2
1.15
1.1

t⋆ increasing

1.05
1
−1

−0.5

0

0.5

1

Z⋆
(a) free swelling

1.45
1.4
1.35

t⋆ increasing

J

1.3
1.25
1.2
1.15
1.1
−1

−0.5

0

0.5

1

Z⋆
(b) compression after free swelling

Figure 5. Ratio of volume of the swollen solid to volume of unswollen solid
(J = det(F s )) as a function of time for (a) free swelling and (b) compressive
force F ⋆ = 1 is applied after the swollen solid reaches a saturated state due
to free swelling. The parameters chosen are χ = 0.425, β1 = 1.3, β2 = 0.018,
µ̄⋆p = 0.1, µ̄⋆G = 0.1, γ ⋆ = 20. Here, 301 spatial points were used for the
calculations, non-dimensional time step chosen is ∆t⋆ = 0.025.
(4) the volume of the mixture is sum of the volumes of the constituents in their natural
state,
(5) the response of the solid from the current configuration to its natural configuration
is isotropic and elastic,
21

1.45
1.4
1.35

J

1.3
1.25
1.2
1.15
1.1
−1

t⋆ increasing

−0.5

0

0.5

1

Z⋆

Figure 6. Ratio of volume of the swollen solid to volume of unswollen solid
(J = det(F s )) as a function of time for free- swelling after the compressive
force is removed. The parameters chosen are χ = 0.425, β1 = 1.3, β2 = 0.018,
µ⋆p = 0.1, µ⋆G = 0.1, γ ⋆ = 20. Here, 301 spatial points were used for the
calculations, non-dimensional time step chosen is ∆t⋆ = 0.025.

(6) the reference configurations of the constituents are same as their natural states,
(7) the rate of dissipation of the mixture is assumed to be the sum of the rates of
dissipation due to mechanical working of the viscoelastic solid, due to the fluid (i.e,
due to the friction between the layers of the fluid), and due to the drag between the
solid and the fluid.
The following additional assumptions are made to solve the problem of compression:
(1) the viscosity of the fluid is zero i.e., we are assuming that the dissipation due to the
friction between the layers of the fluid is much smaller than that due to the drag
between the solid and the fluid,
(2) the velocity and acceleration of the solid are small compared to that of the fluid,
(3) acceleration of the fluid is small,
(4) the external loading is applied on the solid only, whereas the fluid pressure at the
boundary is borne by both the solid and the fluid.

Acknowledgements
The author would like to thank AFOSR/AFRL for their support. Part of this work has
been done when the author was appointed as a lecturer during his Ph.D. by the Department
of Mechanical Engineering at Texas A&M University. This support by the department is
appreciated. The author would also like to thank Mr. Shriram Srinivasan and Dr. K. B.
Nakshatrala of TAMU for their help.
22

References
Atkin, R.J., Craine, R.E., 1976. Continuum theories of mixtures: basic theory and historical
development. The Quarterly Journal of Mechanics and Applied Mathematics 29, 209–244.
Baek, S., Srinivasa, A.R., 2004. Diffusion of a fluid through an elastic solid undergoing large
deformation. International Journal of Non-Linear Mechanics 39, 201–218.
Bedford, A., Drumheller, D.S., 1983. Theories of immiscible and structured mixtures. International Journal of Engineering Science 21, 863–960.
Bhargava, P., 2007. High temperature properties of HFPE-II-52 polyimide resin and composites. Ph.D. thesis. Cornell University.
Biot, M.A., 1956. Theory of deformation of a porous viscoelastic anisotropic solid. Journal
of Applied Physics 27, 459–467.
Bowen, R.M., 1982. Compressible porous media models by use of the theory of mixtures.
International Journal of Engineering Science 20, 697–735.
Cairncross, R.A., Durning, C.J., 1996. A model for drying of viscoelastic polymer coatings.
AIChE Journal 42, 2415–2425.
Cohen, D.S., White Jr., A.B., 1991. Sharp fronts due to diffusion and viscoelastic relaxation
in polymers. SIAM Journal on Applied Mathematics 51, 472–483.
Cohen, S.C., 1992. Postseismic deformation and stress diffusion due to viscoelasticity and
comments on the modified Elsasser model. Journal of Geophysical Research 97, 15395–
15403.
Darcy, H., 1856. Les Fontaines Publiques de la Ville de Dijon. Victor Dalmont, Paris .
Edwards, D.A., Cohen, D.S., 1995. A mathematical model for a dissolving polymer. AIChE
Journal 41, 2345–2355.
Falcone, C.M., Ruggles-Wrenn, M.B., 2009. Rate dependence and short-term creep behavior
of a thermoset polymer at elevated temperature. Journal of Pressure Vessel Technology
131, 011403(1–8).
Fick, A., 1855. Uber diffusion. Poggendorff’s Annalen der Physik und Chemie 94.
Gattiglia, E., Russell, T.P., 1989. Swelling behavior of an aromatic polyimide. Journal of
Polymer Science Part B: Polymer Physics 27, 2131–2144.
Ghosh, M.K., Mittal, K.L. (Eds.), 1996. Polyimides: Fundamentals and Applications. Marcel
Deckker, Inc., New York.
Green, A.E., Naghdi, P.M., 1969. On basic equations for mixtures. The Quarterly Journal
of Mechanics and Applied Mathematics 22, 427–438.
Hron, J., Humphrey, J.D., Rajagopal, K.R., 2002. Material identification of nonlinear solids
infused with a fluid. Mathematics and Mechanics of Solids 7, 629–646.
Huang, S.J., Durning, C.J., 1997. Nonlinear viscoelastic diffusion in concentrated
polystyrene/ethylbenzene solutions. Journal of Polymer Science Part B: Polymer Physics
35, 2103–2119.
Karra, S., Rajagopal, K.R., 2010. Modeling the non-linear viscoelastic response of high
temperature polyimides. Mechanics of Materials In Press, Accepted Manuscript, –.
Kim, Y.R., Little, D.N., Lytton, R.L., 2004. Effect of moisture damage on material properties
and fatigue resistance of asphalt mixtures. Transportation Research Record: Journal of
the Transportation Research Board 1891, 48–54.
Liu, Q., Wang, X., De Kee, D., 2005. Mass transport through swelling membranes. International Journal of Engineering Science 43, 1464–1470.
23

Mills, N., 1966. Incompressible mixtures of Newtonian fluids. International Journal of
Engineering Science 4, 97–112.
Prasad, S.C., Rajagopal, K.R., 2006. On the diffusion of fluids through solids undergoing
large deformations. Mathematics and mechanics of solids 11, 291–305.
Rajagopal, K.R., 1995. Multiple configurations in continuum mechanics. Reports of the
Institute for Computational and Applied Mechanics, University of Pittsburgh (6) .
Rajagopal, K.R., 2003. Diffusion through polymeric solids undergoing large deformations.
Materials Science and Technology 19, 1175–1180.
Rajagopal, K.R., 2007. On a hierarchy of approximate models for flows of incompressible
fluids through porous solids. Mathematical Models and Methods in the Applied Sciences
17, 215–252.
Rajagopal, K.R., Srinivasa, A.R., 2004a. On the thermomechanics of materials that have
multiple natural configurations Part I: Viscoelasticity and classical plasticity. Zeitschrift
für Angewandte Mathematik und Physik (ZAMP) 55, 861–893.
Rajagopal, K.R., Srinivasa, A.R., 2004b. On thermomechanical restrictions of continua. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering
Sciences 460, 631–651.
Rajagopal, K.R., Tao, L., 1995. Mechanics of Mixtures. World Scientific, Farrer Road,
Singapore.
Samohỳl, I., 1987. Thermodynamics of Irreversible Processes in Fluid Mixtures. BG Teubner
Verlagsgesellschaft.
Singh, P.P., Maier, D.E., Cushman, J.H., Campanella, O.H., 2004. Effect of viscoelastic relaxation on moisture transport in foods. Part II: Sorption and drying of soybeans. Journal
of mathematical biology 49, 20–34.
Srinivasan, R., Hall, R.R., Wilson, W.D., Loehle, W.D., Allbee, D.C., 1994. Ultraviolet laser
irradiation of the polyimide, PMDA-ODA (Kapton (TM)), to yield a patternable, porous,
electrically conducting carbon network. Synthetic Metals 66, 301–307.
Truesdell, C., 1957a. Sulle basi della termomeccanica i. Rendiconti Lincei 8, 22–38.
Truesdell, C., 1957b. Sulle basi della termomeccanica ii. Rendiconti Lincei 8, 158–166.
Truesdell, C., Noll, W., Antman, S.S., 2004. The Non-linear Field Theories of Mechanics.
Springer Verlag, Berlin.
Weitsman, Y., 1987. Stress assisted diffusion in elastic and viscoelastic materials. Journal
of the Mechanics and Physics of Solids 35, 73–94.
Yang, C., He, G., He, Y., Ma, P., 2008. Densities and Viscosities of N, N-Dimethylformamide
+ N-Methyl-2-pyrrolidinone and + Dimethyl Sulfoxide in the Temperature Range (303.15
to 353.15) K. Journal of Chemical & Engineering Data 53, 1639–1642.
Appendix A. Convergence of numerical results
Since the analytical solution for the problem is unknown, we perform an engineering
convergence study of the solution using the described algorithm (see figure (7)). In this
study, the error is calculated by taking the difference between solution of various grid sizes
(5, 15, 25, . . . , upto 351 points) and the solution found using a very fine grid of 401 points at
the point Z ⋆ = 0 and at a time of t⋆ = 0.5. Note that the error is propotional to logarithm
of the spatial increment and hence the convergence rate is slow. The aim of our current work
is not to present an optimal algorithm but to solve the coupled partial differential equations.
24

0.02
0.018
0.016

error

0.014
0.012
0.01
0.008
0.006
0.004
0.002

0

1

2

3

4

5

6

-log(∆Z ⋆ )

Figure 7. Engineering spatial convergence of the solution for p at Z ⋆ = 0
and at t⋆ = 0.5. The time-step was chosen to be ∆t⋆ = 0.025.

Appendix B. Derivation of the constitutive equation for the viscoelastic
solid in the absence of diffusion
Here we show that in the absence of diffusion the derived constitutive equations reduces to
a variant of the three-dimensional standard linear solid model given in (Karra and Rajagopal
(2010)). Now, in the absence of diffusion of the fluid, we will have to drop the last term in
(3.33) to get
 
µG0 − µG1 θ
θ
cs1
2
s
s
s
s
+
(IBGs − 3)
ψ̂ = A + (B + c2 ) (θ − θs ) − (θ − θs ) − c2 θln
2
θs
ρs θs
µp0 − µp1 θ
s
+
(IBp(t)
− 3).
(B.1)
ρs θs
We shall assume that the solid is incompressible in the absence of fluid. The constraint of
incompressibility is given by
tr (D s ) = tr (D sG ) = 0.
The reduced energy dissipation equation of the solid reduces to


dψ
s
s
T ·D − ρ
= ξm .
dt θ fixed

(B.2)

(B.3)

In the absence of diffusion, there will be only be dissipation due to mechanical working of
the solid, and so the rate of dissipation in this case would be
ξm = γ(θ)D sG · DsG .
25

(B.4)

Upon maximizing the rate of dissipation using (B.4), (B.2) as constraints (see (Karra and Rajagopal
(2010)) for further details) we arrive at
T s = pI + 2µ̄p B sp(t) ,

(B.5a)

T s = λI + 2µ̄G B sG + ηDsG .

(B.5b)

(p − λ) I + 2µ̄p B sp(t) − 2µ̄G B sG = ηD sG .

(B.6)

From (B.5a), (B.5b) we have
Taking the trace of (B.6), we get

3 (p − λ) = −2µ̄p tr B sp(t) + 2µ̄G tr (B sG ) .

(B.7)

Hence, (B.7) in (B.6) gives



2
µ̄p tr B sp(t) − µ̄G tr (B sG ) + ηDs .
3
The final constitutive equations for the viscoelastic solid are
2µ̄p B sp(t) − 2µ̄G B sG =

T s = pI + 2µ̄p B sp(t) ,

(B.8)
(B.9)

with (B.8) being the evolution equation for the natural configuration of the viscoelastic solid.
This is a variant of the model derived in (Karra and Rajagopal (2010)).
Satish Karra, Texas A&M University, Department of Mechanical Engineering, 3123 TAMU,
College Station TX 77843-3123, United States of America
E-mail address: satkarra@tamu.edu

26

